home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / CMatrix.cc < prev    next >
C/C++ Source or Header  |  1997-07-27  |  80KB  |  3,852 lines
  1. // Matrix manipulations.
  2. /*
  3.  
  4. Copyright (C) 1996 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  21.  
  22. */
  23.  
  24. #if defined (__GNUG__)
  25. #pragma implementation
  26. #endif
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <cfloat>
  33.  
  34. #include <iostream.h>
  35.  
  36. // XXX FIXME XXX
  37. #ifdef HAVE_SYS_TYPES_H
  38. #include <sys/types.h>
  39. #endif
  40.  
  41. #include "CmplxAEPBAL.h"
  42. #include "CmplxDET.h"
  43. #include "CmplxSCHUR.h"
  44. #include "CmplxSVD.h"
  45. #include "f77-fcn.h"
  46. #include "lo-error.h"
  47. #include "lo-ieee.h"
  48. #include "lo-mappers.h"
  49. #include "lo-utils.h"
  50. #include "mx-base.h"
  51. #include "mx-inlines.cc"
  52. #include "oct-cmplx.h"
  53.  
  54. // Fortran functions we call.
  55.  
  56. extern "C"
  57. {
  58.   int F77_FCN (zgemm, ZGEMM) (const char*, const char*, const int&,
  59.                   const int&, const int&, const Complex&,
  60.                   const Complex*, const int&,
  61.                   const Complex*, const int&,
  62.                   const Complex&, Complex*, const int&, 
  63.                   long, long);
  64.  
  65.   int F77_FCN (zgeco, ZGECO) (Complex*, const int&, const int&, int*,
  66.                   double&, Complex*);
  67.  
  68.   int F77_FCN (zgedi, ZGEDI) (Complex*, const int&, const int&, int*,
  69.                   Complex*, Complex*, const int&);
  70.  
  71.   int F77_FCN (zgesl, ZGESL) (Complex*, const int&, const int&, int*,
  72.                   Complex*, const int&);
  73.  
  74.   int F77_FCN (zgelss, ZGELSS) (const int&, const int&, const int&,
  75.                 Complex*, const int&, Complex*,
  76.                 const int&, double*, double&, int&,
  77.                 Complex*, const int&, double*, int&);
  78.  
  79.   // Note that the original complex fft routines were not written for
  80.   // double complex arguments.  They have been modified by adding an
  81.   // implicit double precision (a-h,o-z) statement at the beginning of
  82.   // each subroutine.
  83.  
  84.   int F77_FCN (cffti, CFFTI) (const int&, Complex*);
  85.  
  86.   int F77_FCN (cfftf, CFFTF) (const int&, Complex*, Complex*);
  87.  
  88.   int F77_FCN (cfftb, CFFTB) (const int&, Complex*, Complex*);
  89.  
  90.   int F77_FCN (zlartg, ZLARTG) (const Complex&, const Complex&,
  91.                 double&, Complex&, Complex&);
  92.  
  93.   int F77_FCN (ztrsyl, ZTRSYL) (const char*, const char*, const int&,
  94.                 const int&, const int&,
  95.                 const Complex*, const int&,
  96.                 const Complex*, const int&, 
  97.                 const Complex*, const int&, double&,
  98.                 int&, long, long);
  99.  
  100.   double F77_FCN (zlange, ZLANGE) (const char*, const int&,
  101.                    const int&, const Complex*,
  102.                    const int&, double*); 
  103. }
  104.  
  105. static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
  106.  
  107. // Complex Matrix class
  108.  
  109. ComplexMatrix::ComplexMatrix (const Matrix& a)
  110.   : MArray2<Complex> (a.rows (), a.cols ())
  111. {
  112.   for (int j = 0; j < cols (); j++)
  113.     for (int i = 0; i < rows (); i++)
  114.       elem (i, j) = a.elem (i, j);
  115. }
  116.  
  117. ComplexMatrix::ComplexMatrix (const RowVector& rv)
  118.   : MArray2<Complex> (1, rv.length (), 0.0)
  119. {
  120.   for (int i = 0; i < rv.length (); i++)
  121.     elem (0, i) = rv.elem (i);
  122. }
  123.  
  124. ComplexMatrix::ComplexMatrix (const ColumnVector& cv)
  125.   : MArray2<Complex> (cv.length (), 1, 0.0)
  126. {
  127.   for (int i = 0; i < cv.length (); i++)
  128.     elem (i, 0) = cv.elem (i);
  129. }
  130.  
  131. ComplexMatrix::ComplexMatrix (const DiagMatrix& a)
  132.   : MArray2<Complex> (a.rows (), a.cols (), 0.0)
  133. {
  134.   for (int i = 0; i < a.length (); i++)
  135.     elem (i, i) = a.elem (i, i);
  136. }
  137.  
  138. ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv)
  139.   : MArray2<Complex> (1, rv.length (), 0.0)
  140. {
  141.   for (int i = 0; i < rv.length (); i++)
  142.     elem (0, i) = rv.elem (i);
  143. }
  144.  
  145. ComplexMatrix::ComplexMatrix (const ComplexColumnVector& cv)
  146.   : MArray2<Complex> (cv.length (), 1, 0.0)
  147. {
  148.   for (int i = 0; i < cv.length (); i++)
  149.     elem (i, 0) = cv.elem (i);
  150. }
  151.  
  152. ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a)
  153.   : MArray2<Complex> (a.rows (), a.cols (), 0.0)
  154. {
  155.   for (int i = 0; i < a.length (); i++)
  156.     elem (i, i) = a.elem (i, i);
  157. }
  158.  
  159. // XXX FIXME XXX -- could we use a templated mixed-type copy function
  160. // here?
  161.  
  162. ComplexMatrix::ComplexMatrix (const charMatrix& a)
  163. {
  164.   for (int i = 0; i < a.cols (); i++)
  165.     for (int j = 0; j < a.rows (); j++)
  166.       elem (i, j) = a.elem (i, j);
  167. }
  168.  
  169. bool
  170. ComplexMatrix::operator == (const ComplexMatrix& a) const
  171. {
  172.   if (rows () != a.rows () || cols () != a.cols ())
  173.     return false;
  174.  
  175.   return equal (data (), a.data (), length ());
  176. }
  177.  
  178. bool
  179. ComplexMatrix::operator != (const ComplexMatrix& a) const
  180. {
  181.   return !(*this == a);
  182. }
  183.  
  184. // destructive insert/delete/reorder operations
  185.  
  186. ComplexMatrix&
  187. ComplexMatrix::insert (const Matrix& a, int r, int c)
  188. {
  189.   int a_nr = a.rows ();
  190.   int a_nc = a.cols ();
  191.  
  192.   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
  193.     {
  194.       (*current_liboctave_error_handler) ("range error for insert");
  195.       return *this;
  196.     }
  197.  
  198.   for (int j = 0; j < a_nc; j++)
  199.     for (int i = 0; i < a_nr; i++)
  200.       elem (r+i, c+j) = a.elem (i, j);
  201.  
  202.   return *this;
  203. }
  204.  
  205. ComplexMatrix&
  206. ComplexMatrix::insert (const RowVector& a, int r, int c)
  207. {
  208.   int a_len = a.length ();
  209.   if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
  210.     {
  211.       (*current_liboctave_error_handler) ("range error for insert");
  212.       return *this;
  213.     }
  214.  
  215.   for (int i = 0; i < a_len; i++)
  216.     elem (r, c+i) = a.elem (i);
  217.  
  218.   return *this;
  219. }
  220.  
  221. ComplexMatrix&
  222. ComplexMatrix::insert (const ColumnVector& a, int r, int c)
  223. {
  224.   int a_len = a.length ();
  225.   if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
  226.     {
  227.       (*current_liboctave_error_handler) ("range error for insert");
  228.       return *this;
  229.     }
  230.  
  231.   for (int i = 0; i < a_len; i++)
  232.     elem (r+i, c) = a.elem (i);
  233.  
  234.   return *this;
  235. }
  236.  
  237. ComplexMatrix&
  238. ComplexMatrix::insert (const DiagMatrix& a, int r, int c)
  239. {
  240.   int a_nr = a.rows ();
  241.   int a_nc = a.cols ();
  242.  
  243.   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
  244.     {
  245.       (*current_liboctave_error_handler) ("range error for insert");
  246.       return *this;
  247.     }
  248.  
  249.   fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
  250.  
  251.   for (int i = 0; i < a.length (); i++)
  252.     elem (r+i, c+i) = a.elem (i, i);
  253.  
  254.   return *this;
  255. }
  256.  
  257. ComplexMatrix&
  258. ComplexMatrix::insert (const ComplexMatrix& a, int r, int c)
  259. {
  260.   Array2<Complex>::insert (a, r, c);
  261.   return *this;
  262. }
  263.  
  264. ComplexMatrix&
  265. ComplexMatrix::insert (const ComplexRowVector& a, int r, int c)
  266. {
  267.   int a_len = a.length ();
  268.   if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
  269.     {
  270.       (*current_liboctave_error_handler) ("range error for insert");
  271.       return *this;
  272.     }
  273.  
  274.   for (int i = 0; i < a_len; i++)
  275.     elem (r, c+i) = a.elem (i);
  276.  
  277.   return *this;
  278. }
  279.  
  280. ComplexMatrix&
  281. ComplexMatrix::insert (const ComplexColumnVector& a, int r, int c)
  282. {
  283.   int a_len = a.length ();
  284.   if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
  285.     {
  286.       (*current_liboctave_error_handler) ("range error for insert");
  287.       return *this;
  288.     }
  289.  
  290.   for (int i = 0; i < a_len; i++)
  291.     elem (r+i, c) = a.elem (i);
  292.  
  293.   return *this;
  294. }
  295.  
  296. ComplexMatrix&
  297. ComplexMatrix::insert (const ComplexDiagMatrix& a, int r, int c)
  298. {
  299.   int a_nr = a.rows ();
  300.   int a_nc = a.cols ();
  301.  
  302.   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
  303.     {
  304.       (*current_liboctave_error_handler) ("range error for insert");
  305.       return *this;
  306.     }
  307.  
  308.   fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
  309.  
  310.   for (int i = 0; i < a.length (); i++)
  311.     elem (r+i, c+i) = a.elem (i, i);
  312.  
  313.   return *this;
  314. }
  315.  
  316. ComplexMatrix&
  317. ComplexMatrix::fill (double val)
  318. {
  319.   int nr = rows ();
  320.   int nc = cols ();
  321.   if (nr > 0 && nc > 0)
  322.     for (int j = 0; j < nc; j++)
  323.       for (int i = 0; i < nr; i++)
  324.     elem (i, j) = val;
  325.  
  326.   return *this;
  327. }
  328.  
  329. ComplexMatrix&
  330. ComplexMatrix::fill (const Complex& val)
  331. {
  332.   int nr = rows ();
  333.   int nc = cols ();
  334.   if (nr > 0 && nc > 0)
  335.     for (int j = 0; j < nc; j++)
  336.       for (int i = 0; i < nr; i++)
  337.     elem (i, j) = val;
  338.  
  339.   return *this;
  340. }
  341.  
  342. ComplexMatrix&
  343. ComplexMatrix::fill (double val, int r1, int c1, int r2, int c2)
  344. {
  345.   int nr = rows ();
  346.   int nc = cols ();
  347.   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
  348.       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
  349.     {
  350.       (*current_liboctave_error_handler) ("range error for fill");
  351.       return *this;
  352.     }
  353.  
  354.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  355.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  356.  
  357.   for (int j = c1; j <= c2; j++)
  358.     for (int i = r1; i <= r2; i++)
  359.       elem (i, j) = val;
  360.  
  361.   return *this;
  362. }
  363.  
  364. ComplexMatrix&
  365. ComplexMatrix::fill (const Complex& val, int r1, int c1, int r2, int c2)
  366. {
  367.   int nr = rows ();
  368.   int nc = cols ();
  369.   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
  370.       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
  371.     {
  372.       (*current_liboctave_error_handler) ("range error for fill");
  373.       return *this;
  374.     }
  375.  
  376.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  377.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  378.  
  379.   for (int j = c1; j <= c2; j++)
  380.     for (int i = r1; i <= r2; i++)
  381.       elem (i, j) = val;
  382.  
  383.   return *this;
  384. }
  385.  
  386. ComplexMatrix
  387. ComplexMatrix::append (const Matrix& a) const
  388. {
  389.   int nr = rows ();
  390.   int nc = cols ();
  391.   if (nr != a.rows ())
  392.     {
  393.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  394.       return *this;
  395.     }
  396.  
  397.   int nc_insert = nc;
  398.   ComplexMatrix retval (nr, nc + a.cols ());
  399.   retval.insert (*this, 0, 0);
  400.   retval.insert (a, 0, nc_insert);
  401.   return retval;
  402. }
  403.  
  404. ComplexMatrix
  405. ComplexMatrix::append (const RowVector& a) const
  406. {
  407.   int nr = rows ();
  408.   int nc = cols ();
  409.   if (nr != 1)
  410.     {
  411.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  412.       return *this;
  413.     }
  414.  
  415.   int nc_insert = nc;
  416.   ComplexMatrix retval (nr, nc + a.length ());
  417.   retval.insert (*this, 0, 0);
  418.   retval.insert (a, 0, nc_insert);
  419.   return retval;
  420. }
  421.  
  422. ComplexMatrix
  423. ComplexMatrix::append (const ColumnVector& a) const
  424. {
  425.   int nr = rows ();
  426.   int nc = cols ();
  427.   if (nr != a.length ())
  428.     {
  429.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  430.       return *this;
  431.     }
  432.  
  433.   int nc_insert = nc;
  434.   ComplexMatrix retval (nr, nc + 1);
  435.   retval.insert (*this, 0, 0);
  436.   retval.insert (a, 0, nc_insert);
  437.   return retval;
  438. }
  439.  
  440. ComplexMatrix
  441. ComplexMatrix::append (const DiagMatrix& a) const
  442. {
  443.   int nr = rows ();
  444.   int nc = cols ();
  445.   if (nr != a.rows ())
  446.     {
  447.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  448.       return *this;
  449.     }
  450.  
  451.   int nc_insert = nc;
  452.   ComplexMatrix retval (nr, nc + a.cols ());
  453.   retval.insert (*this, 0, 0);
  454.   retval.insert (a, 0, nc_insert);
  455.   return retval;
  456. }
  457.  
  458. ComplexMatrix
  459. ComplexMatrix::append (const ComplexMatrix& a) const
  460. {
  461.   int nr = rows ();
  462.   int nc = cols ();
  463.   if (nr != a.rows ())
  464.     {
  465.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  466.       return *this;
  467.     }
  468.  
  469.   int nc_insert = nc;
  470.   ComplexMatrix retval (nr, nc + a.cols ());
  471.   retval.insert (*this, 0, 0);
  472.   retval.insert (a, 0, nc_insert);
  473.   return retval;
  474. }
  475.  
  476. ComplexMatrix
  477. ComplexMatrix::append (const ComplexRowVector& a) const
  478. {
  479.   int nr = rows ();
  480.   int nc = cols ();
  481.   if (nr != 1)
  482.     {
  483.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  484.       return *this;
  485.     }
  486.  
  487.   int nc_insert = nc;
  488.   ComplexMatrix retval (nr, nc + a.length ());
  489.   retval.insert (*this, 0, 0);
  490.   retval.insert (a, 0, nc_insert);
  491.   return retval;
  492. }
  493.  
  494. ComplexMatrix
  495. ComplexMatrix::append (const ComplexColumnVector& a) const
  496. {
  497.   int nr = rows ();
  498.   int nc = cols ();
  499.   if (nr != a.length ())
  500.     {
  501.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  502.       return *this;
  503.     }
  504.  
  505.   int nc_insert = nc;
  506.   ComplexMatrix retval (nr, nc + 1);
  507.   retval.insert (*this, 0, 0);
  508.   retval.insert (a, 0, nc_insert);
  509.   return retval;
  510. }
  511.  
  512. ComplexMatrix
  513. ComplexMatrix::append (const ComplexDiagMatrix& a) const
  514. {
  515.   int nr = rows ();
  516.   int nc = cols ();
  517.   if (nr != a.rows ())
  518.     {
  519.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  520.       return *this;
  521.     }
  522.  
  523.   int nc_insert = nc;
  524.   ComplexMatrix retval (nr, nc + a.cols ());
  525.   retval.insert (*this, 0, 0);
  526.   retval.insert (a, 0, nc_insert);
  527.   return retval;
  528. }
  529.  
  530. ComplexMatrix
  531. ComplexMatrix::stack (const Matrix& a) const
  532. {
  533.   int nr = rows ();
  534.   int nc = cols ();
  535.   if (nc != a.cols ())
  536.     {
  537.       (*current_liboctave_error_handler)
  538.     ("column dimension mismatch for stack");
  539.       return *this;
  540.     }
  541.  
  542.   int nr_insert = nr;
  543.   ComplexMatrix retval (nr + a.rows (), nc);
  544.   retval.insert (*this, 0, 0);
  545.   retval.insert (a, nr_insert, 0);
  546.   return retval;
  547. }
  548.  
  549. ComplexMatrix
  550. ComplexMatrix::stack (const RowVector& a) const
  551. {
  552.   int nr = rows ();
  553.   int nc = cols ();
  554.   if (nc != a.length ())
  555.     {
  556.       (*current_liboctave_error_handler)
  557.     ("column dimension mismatch for stack");
  558.       return *this;
  559.     }
  560.  
  561.   int nr_insert = nr;
  562.   ComplexMatrix retval (nr + 1, nc);
  563.   retval.insert (*this, 0, 0);
  564.   retval.insert (a, nr_insert, 0);
  565.   return retval;
  566. }
  567.  
  568. ComplexMatrix
  569. ComplexMatrix::stack (const ColumnVector& a) const
  570. {
  571.   int nr = rows ();
  572.   int nc = cols ();
  573.   if (nc != 1)
  574.     {
  575.       (*current_liboctave_error_handler)
  576.     ("column dimension mismatch for stack");
  577.       return *this;
  578.     }
  579.  
  580.   int nr_insert = nr;
  581.   ComplexMatrix retval (nr + a.length (), nc);
  582.   retval.insert (*this, 0, 0);
  583.   retval.insert (a, nr_insert, 0);
  584.   return retval;
  585. }
  586.  
  587. ComplexMatrix
  588. ComplexMatrix::stack (const DiagMatrix& a) const
  589. {
  590.   int nr = rows ();
  591.   int nc = cols ();
  592.   if (nc != a.cols ())
  593.     {
  594.       (*current_liboctave_error_handler)
  595.     ("column dimension mismatch for stack");
  596.       return *this;
  597.     }
  598.  
  599.   int nr_insert = nr;
  600.   ComplexMatrix retval (nr + a.rows (), nc);
  601.   retval.insert (*this, 0, 0);
  602.   retval.insert (a, nr_insert, 0);
  603.   return retval;
  604. }
  605.  
  606. ComplexMatrix
  607. ComplexMatrix::stack (const ComplexMatrix& a) const
  608. {
  609.   int nr = rows ();
  610.   int nc = cols ();
  611.   if (nc != a.cols ())
  612.     {
  613.       (*current_liboctave_error_handler)
  614.     ("column dimension mismatch for stack");
  615.       return *this;
  616.     }
  617.  
  618.   int nr_insert = nr;
  619.   ComplexMatrix retval (nr + a.rows (), nc);
  620.   retval.insert (*this, 0, 0);
  621.   retval.insert (a, nr_insert, 0);
  622.   return retval;
  623. }
  624.  
  625. ComplexMatrix
  626. ComplexMatrix::stack (const ComplexRowVector& a) const
  627. {
  628.   int nr = rows ();
  629.   int nc = cols ();
  630.   if (nc != a.length ())
  631.     {
  632.       (*current_liboctave_error_handler)
  633.     ("column dimension mismatch for stack");
  634.       return *this;
  635.     }
  636.  
  637.   int nr_insert = nr;
  638.   ComplexMatrix retval (nr + 1, nc);
  639.   retval.insert (*this, 0, 0);
  640.   retval.insert (a, nr_insert, 0);
  641.   return retval;
  642. }
  643.  
  644. ComplexMatrix
  645. ComplexMatrix::stack (const ComplexColumnVector& a) const
  646. {
  647.   int nr = rows ();
  648.   int nc = cols ();
  649.   if (nc != 1)
  650.     {
  651.       (*current_liboctave_error_handler)
  652.     ("column dimension mismatch for stack");
  653.       return *this;
  654.     }
  655.  
  656.   int nr_insert = nr;
  657.   ComplexMatrix retval (nr + a.length (), nc);
  658.   retval.insert (*this, 0, 0);
  659.   retval.insert (a, nr_insert, 0);
  660.   return retval;
  661. }
  662.  
  663. ComplexMatrix
  664. ComplexMatrix::stack (const ComplexDiagMatrix& a) const
  665. {
  666.   int nr = rows ();
  667.   int nc = cols ();
  668.   if (nc != a.cols ())
  669.     {
  670.       (*current_liboctave_error_handler)
  671.     ("column dimension mismatch for stack");
  672.       return *this;
  673.     }
  674.  
  675.   int nr_insert = nr;
  676.   ComplexMatrix retval (nr + a.rows (), nc);
  677.   retval.insert (*this, 0, 0);
  678.   retval.insert (a, nr_insert, 0);
  679.   return retval;
  680. }
  681.  
  682. ComplexMatrix
  683. ComplexMatrix::hermitian (void) const
  684. {
  685.   int nr = rows ();
  686.   int nc = cols ();
  687.   ComplexMatrix result;
  688.   if (length () > 0)
  689.     {
  690.       result.resize (nc, nr);
  691.       for (int j = 0; j < nc; j++)
  692.     for (int i = 0; i < nr; i++)
  693.       result.elem (j, i) = conj (elem (i, j));
  694.     }
  695.   return result;
  696. }
  697.  
  698. ComplexMatrix
  699. ComplexMatrix::transpose (void) const
  700. {
  701.   int nr = rows ();
  702.   int nc = cols ();
  703.   ComplexMatrix result (nc, nr);
  704.   if (length () > 0)
  705.     {
  706.       for (int j = 0; j < nc; j++)
  707.     for (int i = 0; i < nr; i++)
  708.       result.elem (j, i) = elem (i, j);
  709.     }
  710.   return result;
  711. }
  712.  
  713. ComplexMatrix
  714. conj (const ComplexMatrix& a)
  715. {
  716.   int a_len = a.length ();
  717.   ComplexMatrix retval;
  718.   if (a_len > 0)
  719.     retval = ComplexMatrix (conj_dup (a.data (), a_len), a.rows (),
  720.                 a.cols ());
  721.   return retval;
  722. }
  723.  
  724. // resize is the destructive equivalent for this one
  725.  
  726. ComplexMatrix
  727. ComplexMatrix::extract (int r1, int c1, int r2, int c2) const
  728. {
  729.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  730.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  731.  
  732.   int new_r = r2 - r1 + 1;
  733.   int new_c = c2 - c1 + 1;
  734.  
  735.   ComplexMatrix result (new_r, new_c);
  736.  
  737.   for (int j = 0; j < new_c; j++)
  738.     for (int i = 0; i < new_r; i++)
  739.       result.elem (i, j) = elem (r1+i, c1+j);
  740.  
  741.   return result;
  742. }
  743.  
  744. // extract row or column i.
  745.  
  746. ComplexRowVector
  747. ComplexMatrix::row (int i) const
  748. {
  749.   int nc = cols ();
  750.   if (i < 0 || i >= rows ())
  751.     {
  752.       (*current_liboctave_error_handler) ("invalid row selection");
  753.       return ComplexRowVector ();
  754.     }
  755.  
  756.   ComplexRowVector retval (nc);
  757.   for (int j = 0; j < cols (); j++)
  758.     retval.elem (j) = elem (i, j);
  759.  
  760.   return retval;
  761. }
  762.  
  763. ComplexRowVector
  764. ComplexMatrix::row (char *s) const
  765. {
  766.   if (! s)
  767.     {
  768.       (*current_liboctave_error_handler) ("invalid row selection");
  769.       return ComplexRowVector ();
  770.     }
  771.  
  772.   char c = *s;
  773.   if (c == 'f' || c == 'F')
  774.     return row (0);
  775.   else if (c == 'l' || c == 'L')
  776.     return row (rows () - 1);
  777.   else
  778.     {
  779.       (*current_liboctave_error_handler) ("invalid row selection");
  780.       return ComplexRowVector ();
  781.     }
  782. }
  783.  
  784. ComplexColumnVector
  785. ComplexMatrix::column (int i) const
  786. {
  787.   int nr = rows ();
  788.   if (i < 0 || i >= cols ())
  789.     {
  790.       (*current_liboctave_error_handler) ("invalid column selection");
  791.       return ComplexColumnVector ();
  792.     }
  793.  
  794.   ComplexColumnVector retval (nr);
  795.   for (int j = 0; j < nr; j++)
  796.     retval.elem (j) = elem (j, i);
  797.  
  798.   return retval;
  799. }
  800.  
  801. ComplexColumnVector
  802. ComplexMatrix::column (char *s) const
  803. {
  804.   if (! s)
  805.     {
  806.       (*current_liboctave_error_handler) ("invalid column selection");
  807.       return ComplexColumnVector ();
  808.     }
  809.  
  810.   char c = *s;
  811.   if (c == 'f' || c == 'F')
  812.     return column (0);
  813.   else if (c == 'l' || c == 'L')
  814.     return column (cols () - 1);
  815.   else
  816.     {
  817.       (*current_liboctave_error_handler) ("invalid column selection");
  818.       return ComplexColumnVector ();
  819.     }
  820. }
  821.  
  822. ComplexMatrix
  823. ComplexMatrix::inverse (void) const
  824. {
  825.   int info;
  826.   double rcond;
  827.   return inverse (info, rcond);
  828. }
  829.  
  830. ComplexMatrix
  831. ComplexMatrix::inverse (int& info) const
  832. {
  833.   double rcond;
  834.   return inverse (info, rcond);
  835. }
  836.  
  837. ComplexMatrix
  838. ComplexMatrix::inverse (int& info, double& rcond, int force) const
  839. {
  840.   ComplexMatrix retval;
  841.  
  842.   int nr = rows ();
  843.   int nc = cols ();
  844.  
  845.   if (nr != nc)
  846.     (*current_liboctave_error_handler) ("inverse requires square matrix");
  847.   else
  848.     {
  849.       info = 0;
  850.  
  851.       Array<int> ipvt (nr);
  852.       int *pipvt = ipvt.fortran_vec ();
  853.  
  854.       Array<Complex> z (nr);
  855.       Complex *pz = z.fortran_vec ();
  856.  
  857.       retval = *this;
  858.       Complex *tmp_data = retval.fortran_vec ();
  859.  
  860.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nc, pipvt, rcond, pz));
  861.  
  862.       if (f77_exception_encountered)
  863.     (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
  864.       else
  865.     {
  866.       volatile double rcond_plus_one = rcond + 1.0;
  867.  
  868.       if (rcond_plus_one == 1.0)
  869.         info = -1;
  870.  
  871.       if (info == -1 && ! force)
  872.         retval = *this;  // Restore contents.
  873.       else
  874.         {
  875.           Complex *dummy = 0;
  876.  
  877.           F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nc, pipvt, dummy,
  878.                        pz, 1));
  879.  
  880.           if (f77_exception_encountered)
  881.         (*current_liboctave_error_handler)
  882.           ("unrecoverable error in zgedi");
  883.         }
  884.     }
  885.     }
  886.  
  887.   return retval;
  888. }
  889.  
  890. ComplexMatrix
  891. ComplexMatrix::pseudo_inverse (double tol)
  892. {
  893.   ComplexMatrix retval;
  894.  
  895.   ComplexSVD result (*this);
  896.  
  897.   DiagMatrix S = result.singular_values ();
  898.   ComplexMatrix U = result.left_singular_matrix ();
  899.   ComplexMatrix V = result.right_singular_matrix ();
  900.  
  901.   ColumnVector sigma = S.diag ();
  902.  
  903.   int r = sigma.length () - 1;
  904.   int nr = rows ();
  905.   int nc = cols ();
  906.  
  907.   if (tol <= 0.0)
  908.     {
  909.       if (nr > nc)
  910.     tol = nr * sigma.elem (0) * DBL_EPSILON;
  911.       else
  912.     tol = nc * sigma.elem (0) * DBL_EPSILON;
  913.     }
  914.  
  915.   while (r >= 0 && sigma.elem (r) < tol)
  916.     r--;
  917.  
  918.   if (r < 0)
  919.     retval = ComplexMatrix (nc, nr, 0.0);
  920.   else
  921.     {
  922.       ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
  923.       DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
  924.       ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
  925.       retval = Vr * D * Ur.hermitian ();
  926.     }
  927.  
  928.   return retval;
  929. }
  930.  
  931. ComplexMatrix
  932. ComplexMatrix::fourier (void) const
  933. {
  934.   ComplexMatrix retval;
  935.  
  936.   int nr = rows ();
  937.   int nc = cols ();
  938.  
  939.   int npts, nsamples;
  940.  
  941.   if (nr == 1 || nc == 1)
  942.     {
  943.       npts = nr > nc ? nr : nc;
  944.       nsamples = 1;
  945.     }
  946.   else
  947.     {
  948.       npts = nr;
  949.       nsamples = nc;
  950.     }
  951.  
  952.   int nn = 4*npts+15;
  953.  
  954.   Array<Complex> wsave (nn);
  955.   Complex *pwsave = wsave.fortran_vec ();
  956.  
  957.   retval = *this;
  958.   Complex *tmp_data = retval.fortran_vec ();
  959.  
  960.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  961.  
  962.   for (int j = 0; j < nsamples; j++)
  963.     F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
  964.  
  965.   return retval;
  966. }
  967.  
  968. ComplexMatrix
  969. ComplexMatrix::ifourier (void) const
  970. {
  971.   ComplexMatrix retval;
  972.  
  973.   int nr = rows ();
  974.   int nc = cols ();
  975.  
  976.   int npts, nsamples;
  977.  
  978.   if (nr == 1 || nc == 1)
  979.     {
  980.       npts = nr > nc ? nr : nc;
  981.       nsamples = 1;
  982.     }
  983.   else
  984.     {
  985.       npts = nr;
  986.       nsamples = nc;
  987.     }
  988.  
  989.   int nn = 4*npts+15;
  990.  
  991.   Array<Complex> wsave (nn);
  992.   Complex *pwsave = wsave.fortran_vec ();
  993.  
  994.   retval = *this;
  995.   Complex *tmp_data = retval.fortran_vec ();
  996.  
  997.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  998.  
  999.   for (int j = 0; j < nsamples; j++)
  1000.     F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
  1001.  
  1002.   for (int j = 0; j < npts*nsamples; j++)
  1003.     tmp_data[j] = tmp_data[j] / (double) npts;
  1004.  
  1005.   return retval;
  1006. }
  1007.  
  1008. ComplexMatrix
  1009. ComplexMatrix::fourier2d (void) const
  1010. {
  1011.   ComplexMatrix retval;
  1012.  
  1013.   int nr = rows ();
  1014.   int nc = cols ();
  1015.  
  1016.   int npts, nsamples;
  1017.  
  1018.   if (nr == 1 || nc == 1)
  1019.     {
  1020.       npts = nr > nc ? nr : nc;
  1021.       nsamples = 1;
  1022.     }
  1023.   else
  1024.     {
  1025.       npts = nr;
  1026.       nsamples = nc;
  1027.     }
  1028.  
  1029.   int nn = 4*npts+15;
  1030.  
  1031.   Array<Complex> wsave (nn);
  1032.   Complex *pwsave = wsave.fortran_vec ();
  1033.  
  1034.   retval = *this;
  1035.   Complex *tmp_data = retval.fortran_vec ();
  1036.  
  1037.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1038.  
  1039.   for (int j = 0; j < nsamples; j++)
  1040.     F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
  1041.  
  1042.   npts = nc;
  1043.   nsamples = nr;
  1044.   nn = 4*npts+15;
  1045.  
  1046.   wsave.resize (nn);
  1047.   pwsave = wsave.fortran_vec ();
  1048.  
  1049.   Array<Complex> row (npts);
  1050.   Complex *prow = row.fortran_vec ();
  1051.  
  1052.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1053.  
  1054.   for (int j = 0; j < nsamples; j++)
  1055.     {
  1056.       for (int i = 0; i < npts; i++)
  1057.     prow[i] = tmp_data[i*nr + j];
  1058.  
  1059.       F77_FCN (cfftf, CFFTF) (npts, prow, pwsave);
  1060.  
  1061.       for (int i = 0; i < npts; i++)
  1062.     tmp_data[i*nr + j] = prow[i];
  1063.     }
  1064.  
  1065.   return retval;
  1066. }
  1067.  
  1068. ComplexMatrix
  1069. ComplexMatrix::ifourier2d (void) const
  1070. {
  1071.   ComplexMatrix retval;
  1072.  
  1073.   int nr = rows ();
  1074.   int nc = cols ();
  1075.  
  1076.   int npts, nsamples;
  1077.  
  1078.   if (nr == 1 || nc == 1)
  1079.     {
  1080.       npts = nr > nc ? nr : nc;
  1081.       nsamples = 1;
  1082.     }
  1083.   else
  1084.     {
  1085.       npts = nr;
  1086.       nsamples = nc;
  1087.     }
  1088.  
  1089.   int nn = 4*npts+15;
  1090.  
  1091.   Array<Complex> wsave (nn);
  1092.   Complex *pwsave = wsave.fortran_vec ();
  1093.  
  1094.   retval = *this;
  1095.   Complex *tmp_data = retval.fortran_vec ();
  1096.  
  1097.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1098.  
  1099.   for (int j = 0; j < nsamples; j++)
  1100.     F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
  1101.  
  1102.   for (int j = 0; j < npts*nsamples; j++)
  1103.     tmp_data[j] = tmp_data[j] / (double) npts;
  1104.  
  1105.   npts = nc;
  1106.   nsamples = nr;
  1107.   nn = 4*npts+15;
  1108.  
  1109.   wsave.resize (nn);
  1110.   pwsave = wsave.fortran_vec ();
  1111.  
  1112.   Array<Complex> row (npts);
  1113.   Complex *prow = row.fortran_vec ();
  1114.  
  1115.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1116.  
  1117.   for (int j = 0; j < nsamples; j++)
  1118.     {
  1119.       for (int i = 0; i < npts; i++)
  1120.     prow[i] = tmp_data[i*nr + j];
  1121.  
  1122.       F77_FCN (cfftb, CFFTB) (npts, prow, pwsave);
  1123.  
  1124.       for (int i = 0; i < npts; i++)
  1125.     tmp_data[i*nr + j] = prow[i] / (double) npts;
  1126.     }
  1127.  
  1128.   return retval;
  1129. }
  1130.  
  1131. ComplexDET
  1132. ComplexMatrix::determinant (void) const
  1133. {
  1134.   int info;
  1135.   double rcond;
  1136.   return determinant (info, rcond);
  1137. }
  1138.  
  1139. ComplexDET
  1140. ComplexMatrix::determinant (int& info) const
  1141. {
  1142.   double rcond;
  1143.   return determinant (info, rcond);
  1144. }
  1145.  
  1146. ComplexDET
  1147. ComplexMatrix::determinant (int& info, double& rcond) const
  1148. {
  1149.   ComplexDET retval;
  1150.  
  1151.   int nr = rows ();
  1152.   int nc = cols ();
  1153.  
  1154.   if (nr == 0 || nc == 0)
  1155.     {
  1156.       Complex d[2];
  1157.       d[0] = 1.0;
  1158.       d[1] = 0.0;
  1159.       retval = ComplexDET (d);
  1160.     }
  1161.   else
  1162.     {
  1163.       info = 0;
  1164.  
  1165.       Array<int> ipvt (nr);
  1166.       int *pipvt = ipvt.fortran_vec ();
  1167.  
  1168.       Array<Complex> z (nr);
  1169.       Complex *pz = z.fortran_vec ();
  1170.  
  1171.       ComplexMatrix atmp = *this;
  1172.       Complex *tmp_data = atmp.fortran_vec ();
  1173.  
  1174.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  1175.  
  1176.       if (f77_exception_encountered)
  1177.     (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
  1178.       else
  1179.     {
  1180.       volatile double rcond_plus_one = rcond + 1.0;
  1181.  
  1182.       if (rcond_plus_one == 1.0)
  1183.         {
  1184.           info = -1;
  1185.           retval = ComplexDET ();
  1186.         }
  1187.       else
  1188.         {
  1189.           Complex d[2];
  1190.  
  1191.           F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10));
  1192.  
  1193.           if (f77_exception_encountered)
  1194.         (*current_liboctave_error_handler)
  1195.           ("unrecoverable error in dgedi");
  1196.           else
  1197.         retval = ComplexDET (d);
  1198.         }
  1199.     }
  1200.     }
  1201.  
  1202.   return retval;
  1203. }
  1204.  
  1205. ComplexMatrix
  1206. ComplexMatrix::solve (const Matrix& b) const
  1207. {
  1208.   int info;
  1209.   double rcond;
  1210.   return solve (b, info, rcond);
  1211. }
  1212.  
  1213. ComplexMatrix
  1214. ComplexMatrix::solve (const Matrix& b, int& info) const
  1215. {
  1216.   double rcond;
  1217.   return solve (b, info, rcond);
  1218. }
  1219.  
  1220. ComplexMatrix
  1221. ComplexMatrix::solve (const Matrix& b, int& info, double& rcond) const
  1222. {
  1223.   ComplexMatrix tmp (b);
  1224.   return solve (tmp, info, rcond);
  1225. }
  1226.  
  1227. ComplexMatrix
  1228. ComplexMatrix::solve (const ComplexMatrix& b) const
  1229. {
  1230.   int info;
  1231.   double rcond;
  1232.   return solve (b, info, rcond);
  1233. }
  1234.  
  1235. ComplexMatrix
  1236. ComplexMatrix::solve (const ComplexMatrix& b, int& info) const
  1237. {
  1238.   double rcond;
  1239.   return solve (b, info, rcond);
  1240. }
  1241. ComplexMatrix
  1242. ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
  1243. {
  1244.   ComplexMatrix retval;
  1245.  
  1246.   int nr = rows ();
  1247.   int nc = cols ();
  1248.  
  1249.   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
  1250.     (*current_liboctave_error_handler)
  1251.       ("matrix dimension mismatch in solution of linear equations");
  1252.   else
  1253.     {
  1254.       info = 0;
  1255.  
  1256.       Array<int> ipvt (nr);
  1257.       int *pipvt = ipvt.fortran_vec ();
  1258.  
  1259.       Array<Complex> z (nr);
  1260.       Complex *pz = z.fortran_vec ();
  1261.  
  1262.       ComplexMatrix atmp = *this;
  1263.       Complex *tmp_data = atmp.fortran_vec ();
  1264.  
  1265.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  1266.  
  1267.       if (f77_exception_encountered)
  1268.     (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
  1269.       else
  1270.     {
  1271.       volatile double rcond_plus_one = rcond + 1.0;
  1272.  
  1273.       if (rcond_plus_one == 1.0)
  1274.         {
  1275.           info = -2;
  1276.         }
  1277.       else
  1278.         {
  1279.           retval = b;
  1280.           Complex *result = retval.fortran_vec ();
  1281.  
  1282.           int b_nc = b.cols ();
  1283.  
  1284.           for (volatile int j = 0; j < b_nc; j++)
  1285.         {
  1286.           F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt,
  1287.                        &result[nr*j], 0));
  1288.  
  1289.           if (f77_exception_encountered)
  1290.             {
  1291.               (*current_liboctave_error_handler)
  1292.             ("unrecoverable error in dgesl");
  1293.  
  1294.               break;
  1295.             }
  1296.         }
  1297.         }
  1298.     }
  1299.     }
  1300.  
  1301.   return retval;
  1302. }
  1303.  
  1304. ComplexColumnVector
  1305. ComplexMatrix::solve (const ComplexColumnVector& b) const
  1306. {
  1307.   int info;
  1308.   double rcond;
  1309.   return solve (b, info, rcond);
  1310. }
  1311.  
  1312. ComplexColumnVector
  1313. ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const
  1314. {
  1315.   double rcond;
  1316.   return solve (b, info, rcond);
  1317. }
  1318.  
  1319. ComplexColumnVector
  1320. ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
  1321.               double& rcond) const
  1322. {
  1323.   ComplexColumnVector retval;
  1324.  
  1325.   int nr = rows ();
  1326.   int nc = cols ();
  1327.  
  1328.   if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
  1329.     (*current_liboctave_error_handler)
  1330.       ("matrix dimension mismatch in solution of linear equations");
  1331.   else
  1332.     {
  1333.       info = 0;
  1334.  
  1335.       Array<int> ipvt (nr);
  1336.       int *pipvt = ipvt.fortran_vec ();
  1337.  
  1338.       Array<Complex> z (nr);
  1339.       Complex *pz = z.fortran_vec ();
  1340.  
  1341.       ComplexMatrix atmp = *this;
  1342.       Complex *tmp_data = atmp.fortran_vec ();
  1343.  
  1344.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  1345.  
  1346.       if (f77_exception_encountered)
  1347.     (*current_liboctave_error_handler)
  1348.       ("unrecoverable error in dgeco");
  1349.       else
  1350.     {
  1351.       volatile double rcond_plus_one = rcond + 1.0;
  1352.  
  1353.       if (rcond_plus_one == 1.0)
  1354.         {
  1355.           info = -2;
  1356.         }
  1357.       else
  1358.         {
  1359.           retval = b;
  1360.           Complex *result = retval.fortran_vec ();
  1361.  
  1362.           F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, result, 0));
  1363.  
  1364.           if (f77_exception_encountered)
  1365.         (*current_liboctave_error_handler)
  1366.           ("unrecoverable error in dgesl");
  1367.         }
  1368.     }
  1369.     }
  1370.  
  1371.   return retval;
  1372. }
  1373.  
  1374. ComplexMatrix
  1375. ComplexMatrix::lssolve (const ComplexMatrix& b) const
  1376. {
  1377.   int info;
  1378.   int rank;
  1379.   return lssolve (b, info, rank);
  1380. }
  1381.  
  1382. ComplexMatrix
  1383. ComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const
  1384. {
  1385.   int rank;
  1386.   return lssolve (b, info, rank);
  1387. }
  1388.  
  1389. ComplexMatrix
  1390. ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
  1391. {
  1392.   ComplexMatrix retval;
  1393.  
  1394.   int nrhs = b.cols ();
  1395.  
  1396.   int m = rows ();
  1397.   int n = cols ();
  1398.  
  1399.   if (m == 0 || n == 0 || m != b.rows ())
  1400.     (*current_liboctave_error_handler)
  1401.       ("matrix dimension mismatch solution of linear equations");
  1402.   else
  1403.     {
  1404.       ComplexMatrix atmp = *this;
  1405.       Complex *tmp_data = atmp.fortran_vec ();
  1406.  
  1407.       int nrr = m > n ? m : n;
  1408.       ComplexMatrix result (nrr, nrhs);
  1409.  
  1410.       for (int j = 0; j < nrhs; j++)
  1411.     for (int i = 0; i < m; i++)
  1412.       result.elem (i, j) = b.elem (i, j);
  1413.  
  1414.       Complex *presult = result.fortran_vec ();
  1415.  
  1416.       int len_s = m < n ? m : n;
  1417.       Array<double> s (len_s);
  1418.       double *ps = s.fortran_vec ();
  1419.  
  1420.       double rcond = -1.0;
  1421.  
  1422.       int lwork;
  1423.       if (m < n)
  1424.     lwork = 2*m + (nrhs > n ? nrhs : n);
  1425.       else
  1426.     lwork = 2*n + (nrhs > m ? nrhs : m);
  1427.  
  1428.       lwork *= 16;
  1429.  
  1430.       Array<Complex> work (lwork);
  1431.       Complex *pwork = work.fortran_vec ();
  1432.  
  1433.       int lrwork = (5 * (m < n ? m : n)) - 4;
  1434.       lrwork = lrwork > 1 ? lrwork : 1;
  1435.       Array<double> rwork (lrwork);
  1436.       double *prwork = rwork.fortran_vec ();
  1437.  
  1438.       F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
  1439.                  nrr, ps, rcond, rank, pwork, lwork,
  1440.                  prwork, info));
  1441.  
  1442.       if (f77_exception_encountered)
  1443.     (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
  1444.       else
  1445.     {
  1446.       retval.resize (n, nrhs);
  1447.       for (int j = 0; j < nrhs; j++)
  1448.         for (int i = 0; i < n; i++)
  1449.           retval.elem (i, j) = result.elem (i, j);
  1450.     }
  1451.     }
  1452.  
  1453.   return retval;
  1454. }
  1455.  
  1456. ComplexColumnVector
  1457. ComplexMatrix::lssolve (const ComplexColumnVector& b) const
  1458. {
  1459.   int info;
  1460.   int rank;
  1461.   return lssolve (b, info, rank);
  1462. }
  1463.  
  1464. ComplexColumnVector
  1465. ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const
  1466. {
  1467.   int rank;
  1468.   return lssolve (b, info, rank);
  1469. }
  1470.  
  1471. ComplexColumnVector
  1472. ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info,
  1473.             int& rank) const
  1474. {
  1475.   ComplexColumnVector retval;
  1476.  
  1477.   int nrhs = 1;
  1478.  
  1479.   int m = rows ();
  1480.   int n = cols ();
  1481.  
  1482.   if (m == 0 || n == 0 || m != b.length ())
  1483.     (*current_liboctave_error_handler)
  1484.       ("matrix dimension mismatch solution of least squares problem");
  1485.   else
  1486.     {
  1487.       ComplexMatrix atmp = *this;
  1488.       Complex *tmp_data = atmp.fortran_vec ();
  1489.  
  1490.       int nrr = m > n ? m : n;
  1491.       ComplexColumnVector result (nrr);
  1492.  
  1493.       for (int i = 0; i < m; i++)
  1494.     result.elem (i) = b.elem (i);
  1495.  
  1496.       Complex *presult = result.fortran_vec ();
  1497.  
  1498.       int len_s = m < n ? m : n;
  1499.       Array<double> s (len_s);
  1500.       double *ps = s.fortran_vec ();
  1501.  
  1502.       double rcond = -1.0;
  1503.  
  1504.       int lwork;
  1505.       if (m < n)
  1506.     lwork = 2*m + (nrhs > n ? nrhs : n);
  1507.       else
  1508.     lwork = 2*n + (nrhs > m ? nrhs : m);
  1509.  
  1510.       lwork *= 16;
  1511.  
  1512.       Array<Complex> work (lwork);
  1513.       Complex *pwork = work.fortran_vec ();
  1514.  
  1515.       int lrwork = (5 * (m < n ? m : n)) - 4;
  1516.       lrwork = lrwork > 1 ? lrwork : 1;
  1517.       Array<double> rwork (lrwork);
  1518.       double *prwork = rwork.fortran_vec ();
  1519.  
  1520.       F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
  1521.                  nrr, ps, rcond, rank, pwork, lwork,
  1522.                  prwork, info));
  1523.  
  1524.       if (f77_exception_encountered)
  1525.     (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
  1526.       else
  1527.     {
  1528.       retval.resize (n);
  1529.       for (int i = 0; i < n; i++)
  1530.         retval.elem (i) = result.elem (i);
  1531.     }
  1532.     }
  1533.  
  1534.   return retval;
  1535. }
  1536.  
  1537. // Constants for matrix exponential calculation.
  1538.  
  1539. static double padec [] =
  1540. {
  1541.   5.0000000000000000e-1,
  1542.   1.1666666666666667e-1,
  1543.   1.6666666666666667e-2,
  1544.   1.6025641025641026e-3,
  1545.   1.0683760683760684e-4,
  1546.   4.8562548562548563e-6,
  1547.   1.3875013875013875e-7,
  1548.   1.9270852604185938e-9,
  1549. };
  1550.  
  1551. ComplexMatrix
  1552. ComplexMatrix::expm (void) const
  1553. {
  1554.   ComplexMatrix retval;
  1555.  
  1556.   ComplexMatrix m = *this;
  1557.  
  1558.   int nc = columns ();
  1559.  
  1560.   // trace shift value
  1561.   Complex trshift = 0.0;
  1562.  
  1563.   // Preconditioning step 1: trace normalization.
  1564.  
  1565.   for (int i = 0; i < nc; i++)
  1566.     trshift += m.elem (i, i);
  1567.  
  1568.   trshift /= nc;
  1569.  
  1570.   for (int i = 0; i < nc; i++)
  1571.     m.elem (i, i) -= trshift;
  1572.  
  1573.   // Preconditioning step 2: eigenvalue balancing.
  1574.  
  1575.   ComplexAEPBALANCE mbal (m, "B");
  1576.   m = mbal.balanced_matrix ();
  1577.   ComplexMatrix d = mbal.balancing_matrix ();
  1578.  
  1579.   // Preconditioning step 3: scaling.
  1580.  
  1581.   ColumnVector work (nc);
  1582.   double inf_norm
  1583.     = F77_FCN (zlange, ZLANGE) ("I", nc, nc, m.fortran_vec (), nc,
  1584.                 work.fortran_vec ());
  1585.  
  1586.   int sqpow = (int) (inf_norm > 0.0
  1587.              ? (1.0 + log (inf_norm) / log (2.0))
  1588.              : 0.0);
  1589.  
  1590.   // Check whether we need to square at all.
  1591.  
  1592.   if (sqpow < 0)
  1593.     sqpow = 0;
  1594.  
  1595.   if (sqpow > 0)
  1596.     {
  1597.       double scale_factor = 1.0;
  1598.       for (int i = 0; i < sqpow; i++)
  1599.     scale_factor *= 2.0;
  1600.  
  1601.       m = m / scale_factor;
  1602.     }
  1603.  
  1604.   // npp, dpp: pade' approx polynomial matrices.
  1605.  
  1606.   ComplexMatrix npp (nc, nc, 0.0);
  1607.   ComplexMatrix dpp = npp;
  1608.  
  1609.   // Now powers a^8 ... a^1.
  1610.  
  1611.   int minus_one_j = -1;
  1612.   for (int j = 7; j >= 0; j--)
  1613.     {
  1614.       npp = m * npp + m * padec[j];
  1615.       dpp = m * dpp + m * (minus_one_j * padec[j]);
  1616.       minus_one_j *= -1;
  1617.     }
  1618.  
  1619.   // Zero power.
  1620.  
  1621.   dpp = -dpp;
  1622.   for (int j = 0; j < nc; j++)
  1623.     {
  1624.       npp.elem (j, j) += 1.0;
  1625.       dpp.elem (j, j) += 1.0;
  1626.     }
  1627.  
  1628.   // Compute pade approximation = inverse (dpp) * npp.
  1629.  
  1630.   retval = dpp.solve (npp);
  1631.     
  1632.   // Reverse preconditioning step 3: repeated squaring.
  1633.  
  1634.   while (sqpow)
  1635.     {
  1636.       retval = retval * retval;
  1637.       sqpow--;
  1638.     }
  1639.  
  1640.   // Reverse preconditioning step 2: inverse balancing.
  1641.   // XXX FIXME XXX -- should probably do this with Lapack calls
  1642.   // instead of a complete matrix inversion.
  1643.  
  1644.   retval = retval.transpose ();
  1645.   d = d.transpose ();
  1646.   retval = retval * d;
  1647.   retval = d.solve (retval);
  1648.   retval = retval.transpose ();
  1649.  
  1650.   // Reverse preconditioning step 1: fix trace normalization.
  1651.  
  1652.   return retval * exp (trshift);
  1653. }
  1654.  
  1655. // column vector by row vector -> matrix operations
  1656.  
  1657. ComplexMatrix
  1658. operator * (const ColumnVector& v, const ComplexRowVector& a)
  1659. {
  1660.   ComplexColumnVector tmp (v);
  1661.   return tmp * a;
  1662. }
  1663.  
  1664. ComplexMatrix
  1665. operator * (const ComplexColumnVector& a, const RowVector& b)
  1666. {
  1667.   ComplexRowVector tmp (b);
  1668.   return a * tmp;
  1669. }
  1670.  
  1671. ComplexMatrix
  1672. operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
  1673. {
  1674.   ComplexMatrix retval;
  1675.  
  1676.   int len = v.length ();
  1677.   int a_len = a.length ();
  1678.  
  1679.   if (len != a_len)
  1680.     gripe_nonconformant ("operator *", len, 1, 1, a_len);
  1681.   else
  1682.     {
  1683.       if (len != 0)
  1684.     {
  1685.       retval.resize (len, a_len);
  1686.       Complex *c = retval.fortran_vec ();
  1687.  
  1688.       F77_XFCN (zgemm, ZGEMM, ("N", "N", len, a_len, 1, 1.0,
  1689.                    v.data (), len, a.data (), 1, 0.0,
  1690.                    c, len, 1L, 1L)); 
  1691.  
  1692.       if (f77_exception_encountered)
  1693.         (*current_liboctave_error_handler)
  1694.           ("unrecoverable error in zgemm");
  1695.     }
  1696.     }
  1697.  
  1698.   return retval;
  1699. }
  1700.  
  1701. // diagonal matrix by scalar -> matrix operations
  1702.  
  1703. ComplexMatrix
  1704. operator + (const DiagMatrix& a, const Complex& s)
  1705. {
  1706.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1707.   return a + tmp;
  1708. }
  1709.  
  1710. ComplexMatrix
  1711. operator - (const DiagMatrix& a, const Complex& s)
  1712. {
  1713.   ComplexMatrix tmp (a.rows (), a.cols (), -s);
  1714.   return a + tmp;
  1715. }
  1716.  
  1717. ComplexMatrix
  1718. operator + (const ComplexDiagMatrix& a, double s)
  1719. {
  1720.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1721.   return a + tmp;
  1722. }
  1723.  
  1724. ComplexMatrix
  1725. operator - (const ComplexDiagMatrix& a, double s)
  1726. {
  1727.   ComplexMatrix tmp (a.rows (), a.cols (), -s);
  1728.   return a + tmp;
  1729. }
  1730.  
  1731. ComplexMatrix
  1732. operator + (const ComplexDiagMatrix& a, const Complex& s)
  1733. {
  1734.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1735.   return a + tmp;
  1736. }
  1737.  
  1738. ComplexMatrix
  1739. operator - (const ComplexDiagMatrix& a, const Complex& s)
  1740. {
  1741.   ComplexMatrix tmp (a.rows (), a.cols (), -s);
  1742.   return a + tmp;
  1743. }
  1744.  
  1745. // scalar by diagonal matrix -> matrix operations
  1746.  
  1747. ComplexMatrix
  1748. operator + (const Complex& s, const DiagMatrix& a)
  1749. {
  1750.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1751.   return tmp + a;
  1752. }
  1753.  
  1754. ComplexMatrix
  1755. operator - (const Complex& s, const DiagMatrix& a)
  1756. {
  1757.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1758.   return tmp - a;
  1759. }
  1760.  
  1761. ComplexMatrix
  1762. operator + (double s, const ComplexDiagMatrix& a)
  1763. {
  1764.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1765.   return tmp + a;
  1766. }
  1767.  
  1768. ComplexMatrix
  1769. operator - (double s, const ComplexDiagMatrix& a)
  1770. {
  1771.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1772.   return tmp - a;
  1773. }
  1774.  
  1775. ComplexMatrix
  1776. operator + (const Complex& s, const ComplexDiagMatrix& a)
  1777. {
  1778.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1779.   return tmp + a;
  1780. }
  1781.  
  1782. ComplexMatrix
  1783. operator - (const Complex& s, const ComplexDiagMatrix& a)
  1784. {
  1785.   ComplexMatrix tmp (a.rows (), a.cols (), s);
  1786.   return tmp - a;
  1787. }
  1788.  
  1789. // matrix by diagonal matrix -> matrix operations
  1790.  
  1791. ComplexMatrix&
  1792. ComplexMatrix::operator += (const DiagMatrix& a)
  1793. {
  1794.   int nr = rows ();
  1795.   int nc = cols ();
  1796.  
  1797.   int a_nr = rows ();
  1798.   int a_nc = cols ();
  1799.  
  1800.   if (nr != a_nr || nc != a_nc)
  1801.     {
  1802.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1803.       return *this;
  1804.     }
  1805.  
  1806.   for (int i = 0; i < a.length (); i++)
  1807.     elem (i, i) += a.elem (i, i);
  1808.  
  1809.   return *this;
  1810. }
  1811.  
  1812. ComplexMatrix&
  1813. ComplexMatrix::operator -= (const DiagMatrix& a)
  1814. {
  1815.   int nr = rows ();
  1816.   int nc = cols ();
  1817.  
  1818.   int a_nr = rows ();
  1819.   int a_nc = cols ();
  1820.  
  1821.   if (nr != a_nr || nc != a_nc)
  1822.     {
  1823.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1824.       return *this;
  1825.     }
  1826.  
  1827.   for (int i = 0; i < a.length (); i++)
  1828.     elem (i, i) -= a.elem (i, i);
  1829.  
  1830.   return *this;
  1831. }
  1832.  
  1833. ComplexMatrix&
  1834. ComplexMatrix::operator += (const ComplexDiagMatrix& a)
  1835. {
  1836.   int nr = rows ();
  1837.   int nc = cols ();
  1838.  
  1839.   int a_nr = rows ();
  1840.   int a_nc = cols ();
  1841.  
  1842.   if (nr != a_nr || nc != a_nc)
  1843.     {
  1844.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1845.       return *this;
  1846.     }
  1847.  
  1848.   for (int i = 0; i < a.length (); i++)
  1849.     elem (i, i) += a.elem (i, i);
  1850.  
  1851.   return *this;
  1852. }
  1853.  
  1854. ComplexMatrix&
  1855. ComplexMatrix::operator -= (const ComplexDiagMatrix& a)
  1856. {
  1857.   int nr = rows ();
  1858.   int nc = cols ();
  1859.  
  1860.   int a_nr = rows ();
  1861.   int a_nc = cols ();
  1862.  
  1863.   if (nr != a_nr || nc != a_nc)
  1864.     {
  1865.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1866.       return *this;
  1867.     }
  1868.  
  1869.   for (int i = 0; i < a.length (); i++)
  1870.     elem (i, i) -= a.elem (i, i);
  1871.  
  1872.   return *this;
  1873. }
  1874.  
  1875. ComplexMatrix
  1876. operator + (const Matrix& m, const ComplexDiagMatrix& a)
  1877. {
  1878.   int nr = m.rows ();
  1879.   int nc = m.cols ();
  1880.  
  1881.   int a_nr = a.rows ();
  1882.   int a_nc = a.cols ();
  1883.  
  1884.   if (nr != a_nr || nc != a_nc)
  1885.     {
  1886.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  1887.       return ComplexMatrix ();
  1888.     }
  1889.  
  1890.   if (nr == 0 || nc == 0)
  1891.     return ComplexMatrix (nr, nc);
  1892.  
  1893.   ComplexMatrix result (m);
  1894.   for (int i = 0; i < a.length (); i++)
  1895.     result.elem (i, i) += a.elem (i, i);
  1896.  
  1897.   return result;
  1898. }
  1899.  
  1900. ComplexMatrix
  1901. operator - (const Matrix& m, const ComplexDiagMatrix& a)
  1902. {
  1903.   int nr = m.rows ();
  1904.   int nc = m.cols ();
  1905.  
  1906.   int a_nr = a.rows ();
  1907.   int a_nc = a.cols ();
  1908.  
  1909.   if (nr != a_nr || nc != a_nc)
  1910.     {
  1911.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  1912.       return ComplexMatrix ();
  1913.     }
  1914.  
  1915.   if (nr == 0 || nc == 0)
  1916.     return ComplexMatrix (nr, nc);
  1917.  
  1918.   ComplexMatrix result (m);
  1919.   for (int i = 0; i < a.length (); i++)
  1920.     result.elem (i, i) -= a.elem (i, i);
  1921.  
  1922.   return result;
  1923. }
  1924.  
  1925. ComplexMatrix
  1926. operator * (const Matrix& m, const ComplexDiagMatrix& a)
  1927. {
  1928.   ComplexMatrix retval;
  1929.  
  1930.   int nr = m.rows ();
  1931.   int nc = m.cols ();
  1932.  
  1933.   int a_nr = a.rows ();
  1934.   int a_nc = a.cols ();
  1935.  
  1936.   if (nc != a_nr)
  1937.     gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  1938.   else
  1939.     {
  1940.       if (nr == 0 || nc == 0 || a_nc == 0)
  1941.     retval.resize (nr, a_nc, 0.0);
  1942.       else
  1943.     {
  1944.       retval.resize (nr, a_nc);
  1945.       Complex *c = retval.fortran_vec ();
  1946.  
  1947.       Complex *ctmp = 0;
  1948.  
  1949.       for (int j = 0; j < a.length (); j++)
  1950.         {
  1951.           int idx = j * nr;
  1952.           ctmp = c + idx;
  1953.           if (a.elem (j, j) == 1.0)
  1954.         {
  1955.           for (int i = 0; i < nr; i++)
  1956.             ctmp[i] = m.elem (i, j);
  1957.         }
  1958.           else if (a.elem (j, j) == 0.0)
  1959.         {
  1960.           for (int i = 0; i < nr; i++)
  1961.             ctmp[i] = 0.0;
  1962.         }
  1963.           else
  1964.         {
  1965.           for (int i = 0; i < nr; i++)
  1966.             ctmp[i] = a.elem (j, j) * m.elem (i, j);
  1967.         }
  1968.         }
  1969.  
  1970.       if (a_nr < a_nc)
  1971.         {
  1972.           for (int i = nr * nc; i < nr * a_nc; i++)
  1973.         ctmp[i] = 0.0;
  1974.         }
  1975.     }
  1976.     }
  1977.  
  1978.   return retval;
  1979. }
  1980.  
  1981. // diagonal matrix by matrix -> matrix operations
  1982.  
  1983. ComplexMatrix
  1984. operator + (const DiagMatrix& m, const ComplexMatrix& a)
  1985. {
  1986.   int nr = m.rows ();
  1987.   int nc = m.cols ();
  1988.  
  1989.   int a_nr = a.rows ();
  1990.   int a_nc = a.cols ();
  1991.  
  1992.   if (nr != a_nr || nc != a_nc)
  1993.     {
  1994.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  1995.       return ComplexMatrix ();
  1996.     }
  1997.  
  1998.   if (nr == 0 || nc == 0)
  1999.     return ComplexMatrix (nr, nc);
  2000.  
  2001.   ComplexMatrix result (a);
  2002.   for (int i = 0; i < m.length (); i++)
  2003.     result.elem (i, i) += m.elem (i, i);
  2004.  
  2005.   return result;
  2006. }
  2007.  
  2008. ComplexMatrix
  2009. operator - (const DiagMatrix& m, const ComplexMatrix& a)
  2010. {
  2011.   int nr = m.rows ();
  2012.   int nc = m.cols ();
  2013.  
  2014.   int a_nr = a.rows ();
  2015.   int a_nc = a.cols ();
  2016.  
  2017.   if (nr != a_nr || nc != a_nc)
  2018.     {
  2019.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  2020.       return ComplexMatrix ();
  2021.     }
  2022.  
  2023.   if (nr == 0 || nc == 0)
  2024.     return ComplexMatrix (nr, nc);
  2025.  
  2026.   ComplexMatrix result (-a);
  2027.   for (int i = 0; i < m.length (); i++)
  2028.     result.elem (i, i) += m.elem (i, i);
  2029.  
  2030.   return result;
  2031. }
  2032.  
  2033. ComplexMatrix
  2034. operator * (const DiagMatrix& m, const ComplexMatrix& a)
  2035. {
  2036.   int nr = m.rows ();
  2037.   int nc = m.cols ();
  2038.  
  2039.   int a_nr = a.rows ();
  2040.   int a_nc = a.cols ();
  2041.  
  2042.   if (nc != a_nr)
  2043.     {
  2044.       gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  2045.       return ComplexMatrix ();
  2046.     }
  2047.  
  2048.   if (nr == 0 || nc == 0 || a_nc == 0)
  2049.     return ComplexMatrix (nr, nc, 0.0);
  2050.  
  2051.   ComplexMatrix c (nr, a_nc);
  2052.  
  2053.   for (int i = 0; i < m.length (); i++)
  2054.     {
  2055.       if (m.elem (i, i) == 1.0)
  2056.     {
  2057.       for (int j = 0; j < a_nc; j++)
  2058.         c.elem (i, j) = a.elem (i, j);
  2059.     }
  2060.       else if (m.elem (i, i) == 0.0)
  2061.     {
  2062.       for (int j = 0; j < a_nc; j++)
  2063.         c.elem (i, j) = 0.0;
  2064.     }
  2065.       else
  2066.     {
  2067.       for (int j = 0; j < a_nc; j++)
  2068.         c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
  2069.     }
  2070.     }
  2071.  
  2072.   if (nr > nc)
  2073.     {
  2074.       for (int j = 0; j < a_nc; j++)
  2075.     for (int i = a_nr; i < nr; i++)
  2076.       c.elem (i, j) = 0.0;
  2077.     }
  2078.  
  2079.   return c;
  2080. }
  2081.  
  2082. ComplexMatrix
  2083. operator + (const ComplexDiagMatrix& m, const Matrix& a)
  2084. {
  2085.   int nr = m.rows ();
  2086.   int nc = m.cols ();
  2087.  
  2088.   int a_nr = a.rows ();
  2089.   int a_nc = a.cols ();
  2090.  
  2091.   if (nr != a_nr || nc != a_nc)
  2092.     {
  2093.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  2094.       return ComplexMatrix ();
  2095.     }
  2096.  
  2097.   if (nr == 0 || nc == 0)
  2098.     return ComplexMatrix (nr, nc);
  2099.  
  2100.   ComplexMatrix result (a);
  2101.   for (int i = 0; i < m.length (); i++)
  2102.     result.elem (i, i) += m.elem (i, i);
  2103.  
  2104.   return result;
  2105. }
  2106.  
  2107. ComplexMatrix
  2108. operator - (const ComplexDiagMatrix& m, const Matrix& a)
  2109. {
  2110.   int nr = m.rows ();
  2111.   int nc = m.cols ();
  2112.  
  2113.   int a_nr = a.rows ();
  2114.   int a_nc = a.cols ();
  2115.  
  2116.   if (nr != a_nr || nc != a_nc)
  2117.     {
  2118.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  2119.       return ComplexMatrix ();
  2120.     }
  2121.  
  2122.   if (nr == 0 || nc == 0)
  2123.     return ComplexMatrix (nr, nc);
  2124.  
  2125.   ComplexMatrix result (-a);
  2126.   for (int i = 0; i < m.length (); i++)
  2127.     result.elem (i, i) += m.elem (i, i);
  2128.  
  2129.   return result;
  2130. }
  2131.  
  2132. ComplexMatrix
  2133. operator * (const ComplexDiagMatrix& m, const Matrix& a)
  2134. {
  2135.   int nr = m.rows ();
  2136.   int nc = m.cols ();
  2137.  
  2138.   int a_nr = a.rows ();
  2139.   int a_nc = a.cols ();
  2140.  
  2141.   if (nc != a_nr)
  2142.     {
  2143.       gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  2144.       return ComplexMatrix ();
  2145.     }
  2146.  
  2147.   if (nr == 0 || nc == 0 || a_nc == 0)
  2148.     return ComplexMatrix (nr, a_nc, 0.0);
  2149.  
  2150.   ComplexMatrix c (nr, a_nc);
  2151.  
  2152.   for (int i = 0; i < m.length (); i++)
  2153.     {
  2154.       if (m.elem (i, i) == 1.0)
  2155.     {
  2156.       for (int j = 0; j < a_nc; j++)
  2157.         c.elem (i, j) = a.elem (i, j);
  2158.     }
  2159.       else if (m.elem (i, i) == 0.0)
  2160.     {
  2161.       for (int j = 0; j < a_nc; j++)
  2162.         c.elem (i, j) = 0.0;
  2163.     }
  2164.       else
  2165.     {
  2166.       for (int j = 0; j < a_nc; j++)
  2167.         c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
  2168.     }
  2169.     }
  2170.  
  2171.   if (nr > nc)
  2172.     {
  2173.       for (int j = 0; j < a_nc; j++)
  2174.     for (int i = a_nr; i < nr; i++)
  2175.       c.elem (i, j) = 0.0;
  2176.     }
  2177.  
  2178.   return c;
  2179. }
  2180.  
  2181. ComplexMatrix
  2182. operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a)
  2183. {
  2184.   int nr = m.rows ();
  2185.   int nc = m.cols ();
  2186.  
  2187.   int a_nr = a.rows ();
  2188.   int a_nc = a.cols ();
  2189.  
  2190.   if (nr != a_nr || nc != a_nc)
  2191.     {
  2192.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  2193.       return ComplexMatrix ();
  2194.     }
  2195.  
  2196.   if (nr == 0 || nc == 0)
  2197.     return ComplexMatrix (nr, nc);
  2198.  
  2199.   ComplexMatrix result (a);
  2200.   for (int i = 0; i < m.length (); i++)
  2201.     result.elem (i, i) += m.elem (i, i);
  2202.  
  2203.   return result;
  2204. }
  2205.  
  2206. ComplexMatrix
  2207. operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a)
  2208. {
  2209.   int nr = m.rows ();
  2210.   int nc = m.cols ();
  2211.  
  2212.   int a_nr = a.rows ();
  2213.   int a_nc = a.cols ();
  2214.  
  2215.   if (nr != a_nr || nc != a_nc)
  2216.     {
  2217.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  2218.       return ComplexMatrix ();
  2219.     }
  2220.  
  2221.   if (nr == 0 || nc == 0)
  2222.     return ComplexMatrix (nr, nc);
  2223.  
  2224.   ComplexMatrix result (-a);
  2225.   for (int i = 0; i < m.length (); i++)
  2226.     result.elem (i, i) += m.elem (i, i);
  2227.  
  2228.   return result;
  2229. }
  2230.  
  2231. ComplexMatrix
  2232. operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a)
  2233. {
  2234.   int nr = m.rows ();
  2235.   int nc = m.cols ();
  2236.  
  2237.   int a_nr = a.rows ();
  2238.   int a_nc = a.cols ();
  2239.  
  2240.   if (nc != a_nr)
  2241.     {
  2242.       gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  2243.       return ComplexMatrix ();
  2244.     }
  2245.  
  2246.   if (nr == 0 || nc == 0 || a_nc == 0)
  2247.     return ComplexMatrix (nr, a_nc, 0.0);
  2248.  
  2249.   ComplexMatrix c (nr, a_nc);
  2250.  
  2251.   for (int i = 0; i < m.length (); i++)
  2252.     {
  2253.       if (m.elem (i, i) == 1.0)
  2254.     {
  2255.       for (int j = 0; j < a_nc; j++)
  2256.         c.elem (i, j) = a.elem (i, j);
  2257.     }
  2258.       else if (m.elem (i, i) == 0.0)
  2259.     {
  2260.       for (int j = 0; j < a_nc; j++)
  2261.         c.elem (i, j) = 0.0;
  2262.     }
  2263.       else
  2264.     {
  2265.       for (int j = 0; j < a_nc; j++)
  2266.         c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
  2267.     }
  2268.     }
  2269.  
  2270.   if (nr > nc)
  2271.     {
  2272.       for (int j = 0; j < a_nc; j++)
  2273.     for (int i = a_nr; i < nr; i++)
  2274.       c.elem (i, j) = 0.0;
  2275.     }
  2276.  
  2277.   return c;
  2278. }
  2279.  
  2280. // matrix by matrix -> matrix operations
  2281.  
  2282. ComplexMatrix&
  2283. ComplexMatrix::operator += (const Matrix& a)
  2284. {
  2285.   int nr = rows ();
  2286.   int nc = cols ();
  2287.  
  2288.   int a_nr = a.rows ();
  2289.   int a_nc = a.cols ();
  2290.  
  2291.   if (nr != a_nr || nc != a_nc)
  2292.     {
  2293.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  2294.       return *this;
  2295.     }
  2296.  
  2297.   if (nr == 0 || nc == 0)
  2298.     return *this;
  2299.  
  2300.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  2301.  
  2302.   add2 (d, a.data (), length ());
  2303.   return *this;
  2304. }
  2305.  
  2306. ComplexMatrix&
  2307. ComplexMatrix::operator -= (const Matrix& a)
  2308. {
  2309.   int nr = rows ();
  2310.   int nc = cols ();
  2311.  
  2312.   int a_nr = a.rows ();
  2313.   int a_nc = a.cols ();
  2314.  
  2315.   if (nr != a_nr || nc != a_nc)
  2316.     {
  2317.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  2318.       return *this;
  2319.     }
  2320.  
  2321.   if (nr == 0 || nc == 0)
  2322.     return *this;
  2323.  
  2324.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  2325.  
  2326.   subtract2 (d, a.data (), length ());
  2327.   return *this;
  2328. }
  2329.  
  2330. ComplexMatrix&
  2331. ComplexMatrix::operator += (const ComplexMatrix& a)
  2332. {
  2333.   int nr = rows ();
  2334.   int nc = cols ();
  2335.  
  2336.   int a_nr = a.rows ();
  2337.   int a_nc = a.cols ();
  2338.  
  2339.   if (nr != a_nr || nc != a_nc)
  2340.     {
  2341.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  2342.       return *this;
  2343.     }
  2344.  
  2345.   if (nr == 0 || nc == 0)
  2346.     return *this;
  2347.  
  2348.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  2349.  
  2350.   add2 (d, a.data (), length ());
  2351.   return *this;
  2352. }
  2353.  
  2354. ComplexMatrix&
  2355. ComplexMatrix::operator -= (const ComplexMatrix& a)
  2356. {
  2357.   int nr = rows ();
  2358.   int nc = cols ();
  2359.  
  2360.   int a_nr = a.rows ();
  2361.   int a_nc = a.cols ();
  2362.  
  2363.   if (nr != a_nr || nc != a_nc)
  2364.     {
  2365.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  2366.       return *this;
  2367.     }
  2368.  
  2369.   if (nr == 0 || nc == 0)
  2370.     return *this;
  2371.  
  2372.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  2373.  
  2374.   subtract2 (d, a.data (), length ());
  2375.   return *this;
  2376. }
  2377.  
  2378. // unary operations
  2379.  
  2380. Matrix
  2381. ComplexMatrix::operator ! (void) const
  2382. {
  2383.   return Matrix (not (data (), length ()), rows (), cols ());
  2384. }
  2385.  
  2386. // matrix by scalar -> matrix operations
  2387.  
  2388. ComplexMatrix
  2389. operator + (const Matrix& a, const Complex& s)
  2390. {
  2391.   return ComplexMatrix (add (a.data (), a.length (), s),
  2392.             a.rows (), a.cols ());
  2393. }
  2394.  
  2395. ComplexMatrix
  2396. operator - (const Matrix& a, const Complex& s)
  2397. {
  2398.   return ComplexMatrix (subtract (a.data (), a.length (), s),
  2399.             a.rows (), a.cols ());
  2400. }
  2401.  
  2402. ComplexMatrix
  2403. operator * (const Matrix& a, const Complex& s)
  2404. {
  2405.   return ComplexMatrix (multiply (a.data (), a.length (), s),
  2406.             a.rows (), a.cols ());
  2407. }
  2408.  
  2409. ComplexMatrix
  2410. operator / (const Matrix& a, const Complex& s)
  2411. {
  2412.   return ComplexMatrix (divide (a.data (), a.length (), s),
  2413.             a.rows (), a.cols ());
  2414. }
  2415.  
  2416. ComplexMatrix
  2417. operator + (const ComplexMatrix& a, double s)
  2418. {
  2419.   return ComplexMatrix (add (a.data (), a.length (), s),
  2420.             a.rows (), a.cols ());
  2421. }
  2422.  
  2423. ComplexMatrix
  2424. operator - (const ComplexMatrix& a, double s)
  2425. {
  2426.   return ComplexMatrix (subtract (a.data (), a.length (), s),
  2427.             a.rows (), a.cols ());
  2428. }
  2429.  
  2430. ComplexMatrix
  2431. operator * (const ComplexMatrix& a, double s)
  2432. {
  2433.   return ComplexMatrix (multiply (a.data (), a.length (), s),
  2434.             a.rows (), a.cols ());
  2435. }
  2436.  
  2437. ComplexMatrix
  2438. operator / (const ComplexMatrix& a, double s)
  2439. {
  2440.   return ComplexMatrix (divide (a.data (), a.length (), s),
  2441.             a.rows (), a.cols ());
  2442. }
  2443.  
  2444. // scalar by matrix -> matrix operations
  2445.  
  2446. ComplexMatrix
  2447. operator + (double s, const ComplexMatrix& a)
  2448. {
  2449.   return ComplexMatrix (add (a.data (), a.length (), s), a.rows (),
  2450.             a.cols ());
  2451. }
  2452.  
  2453. ComplexMatrix
  2454. operator - (double s, const ComplexMatrix& a)
  2455. {
  2456.   return ComplexMatrix (subtract (s, a.data (), a.length ()),
  2457.             a.rows (), a.cols ());
  2458. }
  2459.  
  2460. ComplexMatrix
  2461. operator * (double s, const ComplexMatrix& a)
  2462. {
  2463.   return ComplexMatrix (multiply (a.data (), a.length (), s),
  2464.             a.rows (), a.cols ());
  2465. }
  2466.  
  2467. ComplexMatrix
  2468. operator / (double s, const ComplexMatrix& a)
  2469. {
  2470.   return ComplexMatrix (divide (s, a.data (), a.length ()),
  2471.             a.rows (), a.cols ());
  2472. }
  2473.  
  2474. ComplexMatrix
  2475. operator + (const Complex& s, const Matrix& a)
  2476. {
  2477.   return ComplexMatrix (add (s, a.data (), a.length ()),
  2478.             a.rows (), a.cols ());
  2479. }
  2480.  
  2481. ComplexMatrix
  2482. operator - (const Complex& s, const Matrix& a)
  2483. {
  2484.   return ComplexMatrix (subtract (s, a.data (), a.length ()),
  2485.             a.rows (), a.cols ());
  2486. }
  2487.  
  2488. ComplexMatrix
  2489. operator * (const Complex& s, const Matrix& a)
  2490. {
  2491.   return ComplexMatrix (multiply (a.data (), a.length (), s),
  2492.             a.rows (), a.cols ());
  2493. }
  2494.  
  2495. ComplexMatrix
  2496. operator / (const Complex& s, const Matrix& a)
  2497. {
  2498.   return ComplexMatrix (divide (s, a.data (), a.length ()),
  2499.             a.rows (), a.cols ());
  2500. }
  2501.  
  2502. // matrix by diagonal matrix -> matrix operations
  2503.  
  2504. ComplexMatrix
  2505. operator + (const ComplexMatrix& m, const DiagMatrix& a)
  2506. {
  2507.   int nr = m.rows ();
  2508.   int nc = m.cols ();
  2509.  
  2510.   int a_nr = a.rows ();
  2511.   int a_nc = a.cols ();
  2512.  
  2513.   if (nr != a_nr || nc != a_nc)
  2514.     {
  2515.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  2516.       return ComplexMatrix ();
  2517.     }
  2518.  
  2519.   if (nr == 0 || nc == 0)
  2520.     return ComplexMatrix (nr, nc);
  2521.  
  2522.   ComplexMatrix result (m);
  2523.   for (int i = 0; i < a.length (); i++)
  2524.     result.elem (i, i) += a.elem (i, i);
  2525.  
  2526.   return result;
  2527. }
  2528.  
  2529. ComplexMatrix
  2530. operator - (const ComplexMatrix& m, const DiagMatrix& a)
  2531. {
  2532.   int nr = m.rows ();
  2533.   int nc = m.cols ();
  2534.  
  2535.   int a_nr = a.rows ();
  2536.   int a_nc = a.cols ();
  2537.  
  2538.   if (nr != a_nr || nc != a_nc)
  2539.     {
  2540.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  2541.       return ComplexMatrix ();
  2542.     }
  2543.  
  2544.   if (nr == 0 || nc == 0)
  2545.     return ComplexMatrix (nr, nc);
  2546.  
  2547.   ComplexMatrix result (m);
  2548.   for (int i = 0; i < a.length (); i++)
  2549.     result.elem (i, i) -= a.elem (i, i);
  2550.  
  2551.   return result;
  2552. }
  2553.  
  2554. ComplexMatrix
  2555. operator * (const ComplexMatrix& m, const DiagMatrix& a)
  2556. {
  2557.   ComplexMatrix retval;
  2558.  
  2559.   int nr = m.rows ();
  2560.   int nc = m.cols ();
  2561.  
  2562.   int a_nr = a.rows ();
  2563.   int a_nc = a.cols ();
  2564.  
  2565.   if (nc != a_nr)
  2566.     gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  2567.   else
  2568.     {
  2569.       if (nr == 0 || nc == 0 || a_nc == 0)
  2570.     retval.resize (nr, nc, 0.0);
  2571.       else
  2572.     {
  2573.       retval.resize (nr, a_nc);
  2574.       Complex *c = retval.fortran_vec ();
  2575.       Complex *ctmp = 0;
  2576.  
  2577.       for (int j = 0; j < a.length (); j++)
  2578.         {
  2579.           int idx = j * nr;
  2580.           ctmp = c + idx;
  2581.           if (a.elem (j, j) == 1.0)
  2582.         {
  2583.           for (int i = 0; i < nr; i++)
  2584.             ctmp[i] = m.elem (i, j);
  2585.         }
  2586.           else if (a.elem (j, j) == 0.0)
  2587.         {
  2588.           for (int i = 0; i < nr; i++)
  2589.             ctmp[i] = 0.0;
  2590.         }
  2591.           else
  2592.         {
  2593.           for (int i = 0; i < nr; i++)
  2594.             ctmp[i] = a.elem (j, j) * m.elem (i, j);
  2595.         }
  2596.         }
  2597.  
  2598.       if (a.rows () < a_nc)
  2599.         {
  2600.           for (int i = nr * nc; i < nr * a_nc; i++)
  2601.         ctmp[i] = 0.0;
  2602.         }
  2603.     }
  2604.     }
  2605.  
  2606.   return retval;
  2607. }
  2608.  
  2609. ComplexMatrix
  2610. operator + (const ComplexMatrix& m, const ComplexDiagMatrix& a)
  2611. {
  2612.   int nr = m.rows ();
  2613.   int nc = m.cols ();
  2614.  
  2615.   int a_nr = a.rows ();
  2616.   int a_nc = a.cols ();
  2617.  
  2618.   if (nr != a_nr || nc != a_nc)
  2619.     {
  2620.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  2621.       return ComplexMatrix ();
  2622.     }
  2623.  
  2624.   if (nr == 0 || nc == 0)
  2625.     return ComplexMatrix (nr, nc);
  2626.  
  2627.   ComplexMatrix result (m);
  2628.   for (int i = 0; i < a.length (); i++)
  2629.     result.elem (i, i) += a.elem (i, i);
  2630.  
  2631.   return result;
  2632. }
  2633.  
  2634. ComplexMatrix
  2635. operator - (const ComplexMatrix& m, const ComplexDiagMatrix& a)
  2636. {
  2637.   int nr = m.rows ();
  2638.   int nc = m.cols ();
  2639.  
  2640.   int a_nr = a.rows ();
  2641.   int a_nc = a.cols ();
  2642.  
  2643.   if (nr != a_nr || nc != a_nc)
  2644.     {
  2645.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  2646.       return ComplexMatrix ();
  2647.     }
  2648.  
  2649.   if (nr == 0 || nc == 0)
  2650.     return ComplexMatrix (nr, nc);
  2651.  
  2652.   ComplexMatrix result (m);
  2653.   for (int i = 0; i < a.length (); i++)
  2654.     result.elem (i, i) -= a.elem (i, i);
  2655.  
  2656.   return result;
  2657. }
  2658.  
  2659. ComplexMatrix
  2660. operator * (const ComplexMatrix& m, const ComplexDiagMatrix& a)
  2661. {
  2662.   ComplexMatrix retval;
  2663.  
  2664.   int nr = m.rows ();
  2665.   int nc = m.cols ();
  2666.  
  2667.   int a_nr = a.rows ();
  2668.   int a_nc = a.cols ();
  2669.  
  2670.   if (nc != a_nr)
  2671.     gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  2672.   else
  2673.     {
  2674.       if (nr == 0 || nc == 0 || a_nc == 0)
  2675.     retval.resize (nr, nc, 0.0);
  2676.       else
  2677.     {
  2678.       retval.resize (nr, nc);
  2679.       Complex *c = retval.fortran_vec ();
  2680.       Complex *ctmp = 0;
  2681.  
  2682.       for (int j = 0; j < a.length (); j++)
  2683.         {
  2684.           int idx = j * nr;
  2685.           ctmp = c + idx;
  2686.           if (a.elem (j, j) == 1.0)
  2687.         {
  2688.           for (int i = 0; i < nr; i++)
  2689.             ctmp[i] = m.elem (i, j);
  2690.         }
  2691.           else if (a.elem (j, j) == 0.0)
  2692.         {
  2693.           for (int i = 0; i < nr; i++)
  2694.             ctmp[i] = 0.0;
  2695.         }
  2696.           else
  2697.         {
  2698.           for (int i = 0; i < nr; i++)
  2699.             ctmp[i] = a.elem (j, j) * m.elem (i, j);
  2700.         }
  2701.         }
  2702.  
  2703.       if (a.rows () < a_nc)
  2704.         {
  2705.           for (int i = nr * nc; i < nr * a_nc; i++)
  2706.         ctmp[i] = 0.0;
  2707.         }
  2708.     }
  2709.     }
  2710.  
  2711.   return retval;
  2712. }
  2713.  
  2714. // matrix by matrix -> matrix operations
  2715.  
  2716. ComplexMatrix
  2717. operator + (const ComplexMatrix& m, const Matrix& a)
  2718. {
  2719.   int nr = m.rows ();
  2720.   int nc = m.cols ();
  2721.  
  2722.   int a_nr = a.rows ();
  2723.   int a_nc = a.cols ();
  2724.  
  2725.   if (nr != a_nr || nc != a_nc)
  2726.     {
  2727.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  2728.       return ComplexMatrix ();
  2729.     }
  2730.  
  2731.   if (nr == 0 || nc == 0)
  2732.     return ComplexMatrix (nr, nc);
  2733.  
  2734.   return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
  2735. }
  2736.  
  2737. ComplexMatrix
  2738. operator - (const ComplexMatrix& m, const Matrix& a)
  2739. {
  2740.   int nr = m.rows ();
  2741.   int nc = m.cols ();
  2742.  
  2743.   int a_nr = a.rows ();
  2744.   int a_nc = a.cols ();
  2745.  
  2746.   if (nr != a_nr || nc != a_nc)
  2747.     {
  2748.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  2749.       return ComplexMatrix ();
  2750.     }
  2751.  
  2752.   if (nr == 0 || nc == 0)
  2753.     return ComplexMatrix (nr, nc);
  2754.  
  2755.   return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
  2756. }
  2757.  
  2758. ComplexMatrix
  2759. operator + (const Matrix& m, const ComplexMatrix& a)
  2760. {
  2761.   int nr = m.rows ();
  2762.   int nc = m.cols ();
  2763.  
  2764.   int a_nr = a.rows ();
  2765.   int a_nc = a.cols ();
  2766.  
  2767.   if (nr != a_nr || nc != a_nc)
  2768.     {
  2769.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  2770.       return ComplexMatrix ();
  2771.     }
  2772.  
  2773.   return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
  2774. }
  2775.  
  2776. ComplexMatrix
  2777. operator - (const Matrix& m, const ComplexMatrix& a)
  2778. {
  2779.   int nr = m.rows ();
  2780.   int nc = m.cols ();
  2781.  
  2782.   int a_nr = a.rows ();
  2783.   int a_nc = a.cols ();
  2784.  
  2785.   if (nr != a_nr || nc != a_nc)
  2786.     {
  2787.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  2788.       return ComplexMatrix ();
  2789.     }
  2790.  
  2791.   if (nr == 0 || nc == 0)
  2792.     return ComplexMatrix (nr, nc);
  2793.  
  2794.   return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
  2795. }
  2796.  
  2797. ComplexMatrix
  2798. operator * (const ComplexMatrix& m, const Matrix& a)
  2799. {
  2800.   ComplexMatrix tmp (a);
  2801.   return m * tmp;
  2802. }
  2803.  
  2804. ComplexMatrix
  2805. operator * (const Matrix& m, const ComplexMatrix& a)
  2806. {
  2807.   ComplexMatrix tmp (m);
  2808.   return tmp * a;
  2809. }
  2810.  
  2811. ComplexMatrix
  2812. operator * (const ComplexMatrix& m, const ComplexMatrix& a)
  2813. {
  2814.   ComplexMatrix retval;
  2815.  
  2816.   int nr = m.rows ();
  2817.   int nc = m.cols ();
  2818.  
  2819.   int a_nr = a.rows ();
  2820.   int a_nc = a.cols ();
  2821.  
  2822.   if (nc != a_nr)
  2823.     gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  2824.   else
  2825.     {
  2826.       if (nr == 0 || nc == 0 || a_nc == 0)
  2827.     retval.resize (nr, nc, 0.0);
  2828.       else
  2829.     {
  2830.       int ld  = nr;
  2831.       int lda = a.rows ();
  2832.  
  2833.       retval.resize (nr, a_nc);
  2834.       Complex *c = retval.fortran_vec ();
  2835.  
  2836.       F77_XFCN (zgemm, ZGEMM, ("N", "N", nr, a_nc, nc, 1.0,
  2837.                    m.data (), ld, a.data (), lda, 0.0,
  2838.                    c, nr, 1L, 1L));
  2839.  
  2840.       if (f77_exception_encountered)
  2841.         (*current_liboctave_error_handler)
  2842.           ("unrecoverable error in zgemm");
  2843.     }
  2844.     }
  2845.  
  2846.   return retval;
  2847. }
  2848.  
  2849. ComplexMatrix
  2850. product (const ComplexMatrix& m, const Matrix& a)
  2851. {
  2852.   int nr = m.rows ();
  2853.   int nc = m.cols ();
  2854.  
  2855.   int a_nr = a.rows ();
  2856.   int a_nc = a.cols ();
  2857.  
  2858.   if (nr != a_nr || nc != a_nc)
  2859.     {
  2860.       gripe_nonconformant ("product", nr, nc, a_nr, a_nc);
  2861.       return ComplexMatrix ();
  2862.     }
  2863.  
  2864.   if (nr == 0 || nc == 0)
  2865.     return ComplexMatrix (nr, nc);
  2866.  
  2867.   return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc);
  2868. }
  2869.  
  2870. ComplexMatrix
  2871. quotient (const ComplexMatrix& m, const Matrix& a)
  2872. {
  2873.   int nr = m.rows ();
  2874.   int nc = m.cols ();
  2875.  
  2876.   int a_nr = a.rows ();
  2877.   int a_nc = a.cols ();
  2878.  
  2879.   if (nr != a_nr || nc != a_nc)
  2880.     {
  2881.       gripe_nonconformant ("quotient", nr, nc, a_nr, a_nc);
  2882.       return ComplexMatrix ();
  2883.     }
  2884.  
  2885.   if (nr == 0 || nc == 0)
  2886.     return ComplexMatrix (nr, nc);
  2887.  
  2888.   return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc);
  2889. }
  2890.  
  2891. ComplexMatrix
  2892. product (const Matrix& m, const ComplexMatrix& a)
  2893. {
  2894.   int nr = m.rows ();
  2895.   int nc = m.cols ();
  2896.  
  2897.   int a_nr = a.rows ();
  2898.   int a_nc = a.cols ();
  2899.  
  2900.   if (nr != a_nr || nc != a_nc)
  2901.     {
  2902.       gripe_nonconformant ("product", nr, nc, a_nr, a_nc);
  2903.       return ComplexMatrix ();
  2904.     }
  2905.  
  2906.   if (nr == 0 || nc == 0)
  2907.     return ComplexMatrix (nr, nc);
  2908.  
  2909.   return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc);
  2910. }
  2911.  
  2912. ComplexMatrix
  2913. quotient (const Matrix& m, const ComplexMatrix& a)
  2914. {
  2915.   int nr = m.rows ();
  2916.   int nc = m.cols ();
  2917.  
  2918.   int a_nr = a.rows ();
  2919.   int a_nc = a.cols ();
  2920.  
  2921.   if (nr != a_nr || nc != a_nc)
  2922.     {
  2923.       gripe_nonconformant ("quotient", nr, nc, a_nr, a_nc);
  2924.       return ComplexMatrix ();
  2925.     }
  2926.  
  2927.   if (nr == 0 || nc == 0)
  2928.     return ComplexMatrix (nr, nc);
  2929.  
  2930.   return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc);
  2931. }
  2932.  
  2933. // other operations
  2934.  
  2935. ComplexMatrix
  2936. ComplexMatrix::map (c_c_Mapper f) const
  2937. {
  2938.   ComplexMatrix b (*this);
  2939.   return b.apply (f);
  2940. }
  2941.  
  2942. Matrix
  2943. ComplexMatrix::map (d_c_Mapper f) const
  2944. {
  2945.   const Complex *d = data ();
  2946.  
  2947.   Matrix retval (rows (), columns ());
  2948.  
  2949.   double *r = retval.fortran_vec ();
  2950.  
  2951.   for (int i = 0; i < length (); i++)
  2952.     r[i] = f (d[i]);
  2953.  
  2954.   return retval;
  2955. }
  2956.  
  2957. ComplexMatrix&
  2958. ComplexMatrix::apply (c_c_Mapper f)
  2959. {
  2960.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  2961.  
  2962.   for (int i = 0; i < length (); i++)
  2963.     d[i] = f (d[i]);
  2964.  
  2965.   return *this;
  2966. }
  2967.  
  2968. bool
  2969. ComplexMatrix::any_element_is_inf_or_nan (void) const
  2970. {
  2971.   int nr = rows ();
  2972.   int nc = cols ();
  2973.  
  2974.   for (int j = 0; j < nc; j++)
  2975.     for (int i = 0; i < nr; i++)
  2976.       {
  2977.     Complex val = elem (i, j);
  2978.     if (xisinf (val) || xisnan (val))
  2979.       return true;
  2980.       }
  2981.  
  2982.   return false;
  2983. }
  2984.  
  2985. // Return true if no elements have imaginary components.
  2986.  
  2987. bool
  2988. ComplexMatrix::all_elements_are_real (void) const
  2989. {
  2990.   int nr = rows ();
  2991.   int nc = cols ();
  2992.  
  2993.   for (int j = 0; j < nc; j++)
  2994.     for (int i = 0; i < nr; i++)
  2995.       if (imag (elem (i, j)) != 0.0)
  2996.     return false;
  2997.  
  2998.   return true;
  2999. }
  3000.  
  3001. // Return nonzero if any element of CM has a non-integer real or
  3002. // imaginary part.  Also extract the largest and smallest (real or
  3003. // imaginary) values and return them in MAX_VAL and MIN_VAL. 
  3004.  
  3005. bool
  3006. ComplexMatrix::all_integers (double& max_val, double& min_val) const
  3007. {
  3008.   int nr = rows ();
  3009.   int nc = cols ();
  3010.  
  3011.   if (nr > 0 && nc > 0)
  3012.     {
  3013.       Complex val = elem (0, 0);
  3014.  
  3015.       double r_val = real (val);
  3016.       double i_val = imag (val);
  3017.  
  3018.       max_val = r_val;
  3019.       min_val = r_val;
  3020.  
  3021.       if (i_val > max_val)
  3022.     max_val = i_val;
  3023.  
  3024.       if (i_val < max_val)
  3025.     min_val = i_val;
  3026.     }
  3027.   else
  3028.     return false;
  3029.  
  3030.   for (int j = 0; j < nc; j++)
  3031.     for (int i = 0; i < nr; i++)
  3032.       {
  3033.     Complex val = elem (i, j);
  3034.  
  3035.     double r_val = real (val);
  3036.     double i_val = imag (val);
  3037.  
  3038.     if (r_val > max_val)
  3039.       max_val = r_val;
  3040.  
  3041.     if (i_val > max_val)
  3042.       max_val = i_val;
  3043.  
  3044.     if (r_val < min_val)
  3045.       min_val = r_val;
  3046.  
  3047.     if (i_val < min_val)
  3048.       min_val = i_val;
  3049.  
  3050.     if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
  3051.       return false;
  3052.       }
  3053.  
  3054.   return true;
  3055. }
  3056.  
  3057. bool
  3058. ComplexMatrix::too_large_for_float (void) const
  3059. {
  3060.   int nr = rows ();
  3061.   int nc = cols ();
  3062.  
  3063.   for (int j = 0; j < nc; j++)
  3064.     for (int i = 0; i < nr; i++)
  3065.       {
  3066.     Complex val = elem (i, j);
  3067.  
  3068.     double r_val = real (val);
  3069.     double i_val = imag (val);
  3070.  
  3071.     if (r_val > FLT_MAX
  3072.         || i_val > FLT_MAX
  3073.         || r_val < FLT_MIN
  3074.         || i_val < FLT_MIN)
  3075.       return true;
  3076.       }
  3077.  
  3078.   return false;
  3079. }
  3080.  
  3081. Matrix
  3082. ComplexMatrix::all (void) const
  3083. {
  3084.   int nr = rows ();
  3085.   int nc = cols ();
  3086.   Matrix retval;
  3087.   if (nr > 0 && nc > 0)
  3088.     {
  3089.       if (nr == 1)
  3090.     {
  3091.       retval.resize (1, 1);
  3092.       retval.elem (0, 0) = 1.0;
  3093.       for (int j = 0; j < nc; j++)
  3094.         {
  3095.           if (elem (0, j) == 0.0)
  3096.         {
  3097.           retval.elem (0, 0) = 0.0;
  3098.           break;
  3099.         }
  3100.         }
  3101.     }
  3102.       else if (nc == 1)
  3103.     {
  3104.       retval.resize (1, 1);
  3105.       retval.elem (0, 0) = 1.0;
  3106.       for (int i = 0; i < nr; i++)
  3107.         {
  3108.           if (elem (i, 0) == 0.0)
  3109.         {
  3110.           retval.elem (0, 0) = 0.0;
  3111.           break;
  3112.         }
  3113.         }
  3114.     }
  3115.       else
  3116.     {
  3117.       retval.resize (1, nc);
  3118.       for (int j = 0; j < nc; j++)
  3119.         {
  3120.           retval.elem (0, j) = 1.0;
  3121.           for (int i = 0; i < nr; i++)
  3122.         {
  3123.           if (elem (i, j) == 0.0)
  3124.             {
  3125.               retval.elem (0, j) = 0.0;
  3126.               break;
  3127.             }
  3128.         }
  3129.         }
  3130.     }
  3131.     }
  3132.   return retval;
  3133. }
  3134.  
  3135. Matrix
  3136. ComplexMatrix::any (void) const
  3137. {
  3138.   int nr = rows ();
  3139.   int nc = cols ();
  3140.   Matrix retval;
  3141.   if (nr > 0 && nc > 0)
  3142.     {
  3143.       if (nr == 1)
  3144.     {
  3145.       retval.resize (1, 1);
  3146.       retval.elem (0, 0) = 0.0;
  3147.       for (int j = 0; j < nc; j++)
  3148.         {
  3149.           if (elem (0, j) != 0.0)
  3150.         {
  3151.           retval.elem (0, 0) = 1.0;
  3152.           break;
  3153.         }
  3154.         }
  3155.     }
  3156.       else if (nc == 1)
  3157.     {
  3158.       retval.resize (1, 1);
  3159.       retval.elem (0, 0) = 0.0;
  3160.       for (int i = 0; i < nr; i++)
  3161.         {
  3162.           if (elem (i, 0) != 0.0)
  3163.         {
  3164.           retval.elem (0, 0) = 1.0;
  3165.           break;
  3166.         }
  3167.         }
  3168.     }
  3169.       else
  3170.     {
  3171.       retval.resize (1, nc);
  3172.       for (int j = 0; j < nc; j++)
  3173.         {
  3174.           retval.elem (0, j) = 0.0;
  3175.           for (int i = 0; i < nr; i++)
  3176.         {
  3177.           if (elem (i, j) != 0.0)
  3178.             {
  3179.               retval.elem (0, j) = 1.0;
  3180.               break;
  3181.             }
  3182.         }
  3183.         }
  3184.     }
  3185.     }
  3186.   return retval;
  3187. }
  3188.  
  3189. ComplexMatrix
  3190. ComplexMatrix::cumprod (void) const
  3191. {
  3192.   int nr = rows ();
  3193.   int nc = cols ();
  3194.   ComplexMatrix retval;
  3195.   if (nr > 0 && nc > 0)
  3196.     {
  3197.       if (nr == 1)
  3198.     {
  3199.       retval.resize (1, nc);
  3200.       Complex prod = elem (0, 0);
  3201.       for (int j = 0; j < nc; j++)
  3202.         {
  3203.           retval.elem (0, j) = prod;
  3204.           if (j < nc - 1)
  3205.         prod *= elem (0, j+1);
  3206.         }
  3207.     }
  3208.       else if (nc == 1)
  3209.     {
  3210.       retval.resize (nr, 1);
  3211.       Complex prod = elem (0, 0);
  3212.       for (int i = 0; i < nr; i++)
  3213.         {
  3214.           retval.elem (i, 0) = prod;
  3215.           if (i < nr - 1)
  3216.         prod *= elem (i+1, 0);
  3217.         }
  3218.     }
  3219.       else
  3220.     {
  3221.       retval.resize (nr, nc);
  3222.       for (int j = 0; j < nc; j++)
  3223.         {
  3224.           Complex prod = elem (0, j);
  3225.           for (int i = 0; i < nr; i++)
  3226.         {
  3227.           retval.elem (i, j) = prod;
  3228.           if (i < nr - 1)
  3229.             prod *= elem (i+1, j);
  3230.         }
  3231.         }
  3232.     }
  3233.     }
  3234.   return retval;
  3235. }
  3236.  
  3237. ComplexMatrix
  3238. ComplexMatrix::cumsum (void) const
  3239. {
  3240.   int nr = rows ();
  3241.   int nc = cols ();
  3242.   ComplexMatrix retval;
  3243.   if (nr > 0 && nc > 0)
  3244.     {
  3245.       if (nr == 1)
  3246.     {
  3247.       retval.resize (1, nc);
  3248.       Complex sum = elem (0, 0);
  3249.       for (int j = 0; j < nc; j++)
  3250.         {
  3251.           retval.elem (0, j) = sum;
  3252.           if (j < nc - 1)
  3253.         sum += elem (0, j+1);
  3254.         }
  3255.     }
  3256.       else if (nc == 1)
  3257.     {
  3258.       retval.resize (nr, 1);
  3259.       Complex sum = elem (0, 0);
  3260.       for (int i = 0; i < nr; i++)
  3261.         {
  3262.           retval.elem (i, 0) = sum;
  3263.           if (i < nr - 1)
  3264.         sum += elem (i+1, 0);
  3265.         }
  3266.     }
  3267.       else
  3268.     {
  3269.       retval.resize (nr, nc);
  3270.       for (int j = 0; j < nc; j++)
  3271.         {
  3272.           Complex sum = elem (0, j);
  3273.           for (int i = 0; i < nr; i++)
  3274.         {
  3275.           retval.elem (i, j) = sum;
  3276.           if (i < nr - 1)
  3277.             sum += elem (i+1, j);
  3278.         }
  3279.         }
  3280.     }
  3281.     }
  3282.   return retval;
  3283. }
  3284.  
  3285. ComplexMatrix
  3286. ComplexMatrix::prod (void) const
  3287. {
  3288.   int nr = rows ();
  3289.   int nc = cols ();
  3290.   ComplexMatrix retval;
  3291.   if (nr > 0 && nc > 0)
  3292.     {
  3293.       if (nr == 1)
  3294.     {
  3295.       retval.resize (1, 1);
  3296.       retval.elem (0, 0) = 1.0;
  3297.       for (int j = 0; j < nc; j++)
  3298.         retval.elem (0, 0) *= elem (0, j);
  3299.     }
  3300.       else if (nc == 1)
  3301.     {
  3302.       retval.resize (1, 1);
  3303.       retval.elem (0, 0) = 1.0;
  3304.       for (int i = 0; i < nr; i++)
  3305.         retval.elem (0, 0) *= elem (i, 0);
  3306.     }
  3307.       else
  3308.     {
  3309.       retval.resize (1, nc);
  3310.       for (int j = 0; j < nc; j++)
  3311.         {
  3312.           retval.elem (0, j) = 1.0;
  3313.           for (int i = 0; i < nr; i++)
  3314.         retval.elem (0, j) *= elem (i, j);
  3315.         }
  3316.     }
  3317.     }
  3318.   return retval;
  3319. }
  3320.  
  3321. ComplexMatrix
  3322. ComplexMatrix::sum (void) const
  3323. {
  3324.   int nr = rows ();
  3325.   int nc = cols ();
  3326.   ComplexMatrix retval;
  3327.   if (nr > 0 && nc > 0)
  3328.     {
  3329.       if (nr == 1)
  3330.     {
  3331.       retval.resize (1, 1);
  3332.       retval.elem (0, 0) = 0.0;
  3333.       for (int j = 0; j < nc; j++)
  3334.         retval.elem (0, 0) += elem (0, j);
  3335.     }
  3336.       else if (nc == 1)
  3337.     {
  3338.       retval.resize (1, 1);
  3339.       retval.elem (0, 0) = 0.0;
  3340.       for (int i = 0; i < nr; i++)
  3341.         retval.elem (0, 0) += elem (i, 0);
  3342.     }
  3343.       else
  3344.     {
  3345.       retval.resize (1, nc);
  3346.       for (int j = 0; j < nc; j++)
  3347.         {
  3348.           retval.elem (0, j) = 0.0;
  3349.           for (int i = 0; i < nr; i++)
  3350.         retval.elem (0, j) += elem (i, j);
  3351.         }
  3352.     }
  3353.     }
  3354.   return retval;
  3355. }
  3356.  
  3357. ComplexMatrix
  3358. ComplexMatrix::sumsq (void) const
  3359. {
  3360.   int nr = rows ();
  3361.   int nc = cols ();
  3362.   ComplexMatrix retval;
  3363.   if (nr > 0 && nc > 0)
  3364.     {
  3365.       if (nr == 1)
  3366.     {
  3367.       retval.resize (1, 1);
  3368.       retval.elem (0, 0) = 0.0;
  3369.       for (int j = 0; j < nc; j++)
  3370.         {
  3371.           Complex d = elem (0, j);
  3372.           retval.elem (0, 0) += d * d;
  3373.         }
  3374.     }
  3375.       else if (nc == 1)
  3376.     {
  3377.       retval.resize (1, 1);
  3378.       retval.elem (0, 0) = 0.0;
  3379.       for (int i = 0; i < nr; i++)
  3380.         {
  3381.           Complex d = elem (i, 0);
  3382.           retval.elem (0, 0) += d * d;
  3383.         }
  3384.     }
  3385.       else
  3386.     {
  3387.       retval.resize (1, nc);
  3388.       for (int j = 0; j < nc; j++)
  3389.         {
  3390.           retval.elem (0, j) = 0.0;
  3391.           for (int i = 0; i < nr; i++)
  3392.         {
  3393.           Complex d = elem (i, j);
  3394.           retval.elem (0, j) += d * d;
  3395.         }
  3396.         }
  3397.     }
  3398.     }
  3399.   return retval;
  3400. }
  3401.  
  3402. ComplexColumnVector
  3403. ComplexMatrix::diag (void) const
  3404. {
  3405.   return diag (0);
  3406. }
  3407.  
  3408. ComplexColumnVector
  3409. ComplexMatrix::diag (int k) const
  3410. {
  3411.   int nnr = rows ();
  3412.   int nnc = cols ();
  3413.   if (k > 0)
  3414.     nnc -= k;
  3415.   else if (k < 0)
  3416.     nnr += k;
  3417.  
  3418.   ComplexColumnVector d;
  3419.  
  3420.   if (nnr > 0 && nnc > 0)
  3421.     {
  3422.       int ndiag = (nnr < nnc) ? nnr : nnc;
  3423.  
  3424.       d.resize (ndiag);
  3425.  
  3426.       if (k > 0)
  3427.     {
  3428.       for (int i = 0; i < ndiag; i++)
  3429.         d.elem (i) = elem (i, i+k);
  3430.     }
  3431.       else if ( k < 0)
  3432.     {
  3433.       for (int i = 0; i < ndiag; i++)
  3434.         d.elem (i) = elem (i-k, i);
  3435.     }
  3436.       else
  3437.     {
  3438.       for (int i = 0; i < ndiag; i++)
  3439.         d.elem (i) = elem (i, i);
  3440.     }
  3441.     }
  3442.   else
  3443.     cerr << "diag: requested diagonal out of range\n";
  3444.  
  3445.   return d;
  3446. }
  3447.  
  3448. bool
  3449. ComplexMatrix::row_is_real_only (int i) const
  3450. {
  3451.   bool retval = true;
  3452.  
  3453.   int nc = columns ();
  3454.  
  3455.   for (int j = 0; j < nc; j++)
  3456.     {
  3457.       if (imag (elem (i, j)) != 0.0)
  3458.     {
  3459.       retval = false;
  3460.       break;
  3461.     }
  3462.     }
  3463.  
  3464.   return retval;          
  3465. }
  3466.  
  3467. bool
  3468. ComplexMatrix::column_is_real_only (int j) const
  3469. {
  3470.   bool retval = true;
  3471.  
  3472.   int nr = rows ();
  3473.  
  3474.   for (int i = 0; i < nr; i++)
  3475.     {
  3476.       if (imag (elem (i, j)) != 0.0)
  3477.     {
  3478.       retval = false;
  3479.       break;
  3480.     }
  3481.     }
  3482.  
  3483.   return retval;          
  3484. }
  3485.  
  3486. ComplexColumnVector
  3487. ComplexMatrix::row_min (void) const
  3488. {
  3489.   Array<int> index;
  3490.   return row_min (index);
  3491. }
  3492.  
  3493. ComplexColumnVector
  3494. ComplexMatrix::row_min (Array<int>& index) const
  3495. {
  3496.   ComplexColumnVector result;
  3497.  
  3498.   int nr = rows ();
  3499.   int nc = cols ();
  3500.  
  3501.   if (nr > 0 && nc > 0)
  3502.     {
  3503.       result.resize (nr);
  3504.       index.resize (nr);
  3505.  
  3506.       for (int i = 0; i < nr; i++)
  3507.         {
  3508.       int idx = 0;
  3509.  
  3510.       Complex tmp_min = elem (i, idx);
  3511.  
  3512.       bool real_only = row_is_real_only (i);
  3513.  
  3514.       double abs_min = real_only ? real (tmp_min) : abs (tmp_min);
  3515.  
  3516.       if (xisnan (tmp_min))
  3517.         idx = -1;
  3518.       else
  3519.         {
  3520.           for (int j = 1; j < nc; j++)
  3521.         {
  3522.           Complex tmp = elem (i, j);
  3523.  
  3524.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  3525.  
  3526.           if (xisnan (tmp))
  3527.             {
  3528.               idx = -1;
  3529.               break;
  3530.             }
  3531.           else if (abs_tmp < abs_min)
  3532.             {
  3533.               idx = j;
  3534.               tmp_min = tmp;
  3535.               abs_min = abs_tmp;
  3536.             }
  3537.         }
  3538.  
  3539.           result.elem (i) = (idx < 0) ? Complex_NaN_result : tmp_min;
  3540.           index.elem (i) = idx;
  3541.         }
  3542.         }
  3543.     }
  3544.  
  3545.   return result;
  3546. }
  3547.  
  3548. ComplexColumnVector
  3549. ComplexMatrix::row_max (void) const
  3550. {
  3551.   Array<int> index;
  3552.   return row_max (index);
  3553. }
  3554.  
  3555. ComplexColumnVector
  3556. ComplexMatrix::row_max (Array<int>& index) const
  3557. {
  3558.   ComplexColumnVector result;
  3559.  
  3560.   int nr = rows ();
  3561.   int nc = cols ();
  3562.  
  3563.   if (nr > 0 && nc > 0)
  3564.     {
  3565.       result.resize (nr);
  3566.       index.resize (nr);
  3567.  
  3568.       for (int i = 0; i < nr; i++)
  3569.         {
  3570.       int idx = 0;
  3571.  
  3572.       Complex tmp_max = elem (i, idx);
  3573.  
  3574.       bool real_only = row_is_real_only (i);
  3575.  
  3576.       double abs_max = real_only ? real (tmp_max) : abs (tmp_max);
  3577.  
  3578.       if (xisnan (tmp_max))
  3579.         idx = -1;
  3580.       else
  3581.         {
  3582.           for (int j = 1; j < nc; j++)
  3583.         {
  3584.           Complex tmp = elem (i, j);
  3585.  
  3586.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  3587.  
  3588.           if (xisnan (tmp))
  3589.             {
  3590.               idx = -1;
  3591.               break;
  3592.             }
  3593.           else if (abs_tmp > abs_max)
  3594.             {
  3595.               idx = j;
  3596.               tmp_max = tmp;
  3597.               abs_max = abs_tmp;
  3598.             }
  3599.         }
  3600.  
  3601.           result.elem (i) = (idx < 0) ? Complex_NaN_result : tmp_max;
  3602.           index.elem (i) = idx;
  3603.         }
  3604.         }
  3605.     }
  3606.  
  3607.   return result;
  3608. }
  3609.  
  3610. ComplexRowVector
  3611. ComplexMatrix::column_min (void) const
  3612. {
  3613.   Array<int> index;
  3614.   return column_min (index);
  3615. }
  3616.  
  3617. ComplexRowVector
  3618. ComplexMatrix::column_min (Array<int>& index) const
  3619. {
  3620.   ComplexRowVector result;
  3621.  
  3622.   int nr = rows ();
  3623.   int nc = cols ();
  3624.  
  3625.   if (nr > 0 && nc > 0)
  3626.     {
  3627.       result.resize (nc);
  3628.       index.resize (nc);
  3629.  
  3630.       for (int j = 0; j < nc; j++)
  3631.         {
  3632.       int idx = 0;
  3633.  
  3634.       Complex tmp_min = elem (idx, j);
  3635.  
  3636.       bool real_only = column_is_real_only (j);
  3637.  
  3638.       double abs_min = real_only ? real (tmp_min) : abs (tmp_min);
  3639.  
  3640.       if (xisnan (tmp_min))
  3641.         idx = -1;
  3642.       else
  3643.         {
  3644.           for (int i = 1; i < nr; i++)
  3645.         {
  3646.           Complex tmp = elem (i, j);
  3647.  
  3648.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  3649.  
  3650.           if (xisnan (tmp))
  3651.             {
  3652.               idx = -1;
  3653.               break;
  3654.             }
  3655.           else if (abs_tmp < abs_min)
  3656.             {
  3657.               idx = i;
  3658.               tmp_min = tmp;
  3659.               abs_min = abs_tmp;
  3660.             }
  3661.         }
  3662.  
  3663.           result.elem (j) = (idx < 0) ? Complex_NaN_result : tmp_min;
  3664.           index.elem (j) = idx;
  3665.         }
  3666.         }
  3667.     }
  3668.  
  3669.   return result;
  3670. }
  3671.  
  3672. ComplexRowVector
  3673. ComplexMatrix::column_max (void) const
  3674. {
  3675.   Array<int> index;
  3676.   return column_max (index);
  3677. }
  3678.  
  3679. ComplexRowVector
  3680. ComplexMatrix::column_max (Array<int>& index) const
  3681. {
  3682.   ComplexRowVector result;
  3683.  
  3684.   int nr = rows ();
  3685.   int nc = cols ();
  3686.  
  3687.   if (nr > 0 && nc > 0)
  3688.     {
  3689.       result.resize (nc);
  3690.       index.resize (nc);
  3691.  
  3692.       for (int j = 0; j < nc; j++)
  3693.         {
  3694.       int idx = 0;
  3695.  
  3696.       Complex tmp_max = elem (idx, j);
  3697.  
  3698.       bool real_only = column_is_real_only (j);
  3699.  
  3700.       double abs_max = real_only ? real (tmp_max) : abs (tmp_max);
  3701.  
  3702.       if (xisnan (tmp_max))
  3703.         idx = -1;
  3704.       else
  3705.         {
  3706.           for (int i = 1; i < nr; i++)
  3707.         {
  3708.           Complex tmp = elem (i, j);
  3709.  
  3710.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  3711.  
  3712.           if (xisnan (tmp))
  3713.             {
  3714.               idx = -1;
  3715.               break;
  3716.             }
  3717.           else if (abs_tmp > abs_max)
  3718.             {
  3719.               idx = i;
  3720.               tmp_max = tmp;
  3721.               abs_max = abs_tmp;
  3722.             }
  3723.         }
  3724.  
  3725.           result.elem (j) = (idx < 0) ? Complex_NaN_result : tmp_max;
  3726.           index.elem (j) = idx;
  3727.         }
  3728.         }
  3729.     }
  3730.  
  3731.   return result;
  3732. }
  3733.  
  3734. // i/o
  3735.  
  3736. ostream&
  3737. operator << (ostream& os, const ComplexMatrix& a)
  3738. {
  3739. //  int field_width = os.precision () + 7;
  3740.   for (int i = 0; i < a.rows (); i++)
  3741.     {
  3742.       for (int j = 0; j < a.cols (); j++)
  3743.     os << " " /* setw (field_width) */ << a.elem (i, j);
  3744.       os << "\n";
  3745.     }
  3746.   return os;
  3747. }
  3748.  
  3749. istream&
  3750. operator >> (istream& is, ComplexMatrix& a)
  3751. {
  3752.   int nr = a.rows ();
  3753.   int nc = a.cols ();
  3754.  
  3755.   if (nr < 1 || nc < 1)
  3756.     is.clear (ios::badbit);
  3757.   else
  3758.     {
  3759.       Complex tmp;
  3760.       for (int i = 0; i < nr; i++)
  3761.     for (int j = 0; j < nc; j++)
  3762.       {
  3763.         is >> tmp;
  3764.         if (is)
  3765.           a.elem (i, j) = tmp;
  3766.         else
  3767.           goto done;
  3768.       }
  3769.     }
  3770.  
  3771. done:
  3772.  
  3773.   return is;
  3774. }
  3775.  
  3776. ComplexMatrix
  3777. Givens (const Complex& x, const Complex& y)
  3778. {
  3779.   double cc;
  3780.   Complex cs, temp_r;
  3781.  
  3782.   F77_FCN (zlartg, ZLARTG) (x, y, cc, cs, temp_r);
  3783.  
  3784.   ComplexMatrix g (2, 2);
  3785.  
  3786.   g.elem (0, 0) = cc;
  3787.   g.elem (1, 1) = cc;
  3788.   g.elem (0, 1) = cs;
  3789.   g.elem (1, 0) = -conj (cs);
  3790.  
  3791.   return g;
  3792. }
  3793.  
  3794. ComplexMatrix
  3795. Sylvester (const ComplexMatrix& a, const ComplexMatrix& b,
  3796.        const ComplexMatrix& c)
  3797. {
  3798.   ComplexMatrix retval;
  3799.  
  3800.   // XXX FIXME XXX -- need to check that a, b, and c are all the same
  3801.   // size.
  3802.  
  3803.   // Compute Schur decompositions
  3804.  
  3805.   ComplexSCHUR as (a, "U");
  3806.   ComplexSCHUR bs (b, "U");
  3807.   
  3808.   // Transform c to new coordinates.
  3809.  
  3810.   ComplexMatrix ua = as.unitary_matrix ();
  3811.   ComplexMatrix sch_a = as.schur_matrix ();
  3812.  
  3813.   ComplexMatrix ub = bs.unitary_matrix ();
  3814.   ComplexMatrix sch_b = bs.schur_matrix ();
  3815.   
  3816.   ComplexMatrix cx = ua.hermitian () * c * ub;
  3817.  
  3818.   // Solve the sylvester equation, back-transform, and return the
  3819.   // solution.
  3820.  
  3821.   int a_nr = a.rows ();
  3822.   int b_nr = b.rows ();
  3823.  
  3824.   double scale;
  3825.   int info;
  3826.  
  3827.   Complex *pa = sch_a.fortran_vec ();
  3828.   Complex *pb = sch_b.fortran_vec ();
  3829.   Complex *px = cx.fortran_vec ();
  3830.   
  3831.   F77_XFCN (ztrsyl, ZTRSYL, ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb,
  3832.                  b_nr, px, a_nr, scale,
  3833.                  info, 1L, 1L));
  3834.  
  3835.   if (f77_exception_encountered)
  3836.     (*current_liboctave_error_handler) ("unrecoverable error in ztrsyl");
  3837.   else
  3838.     {
  3839.       // XXX FIXME XXX -- check info?
  3840.  
  3841.       retval = -ua * cx * ub.hermitian ();
  3842.     }
  3843.  
  3844.   return retval;
  3845. }
  3846.  
  3847. /*
  3848. ;;; Local Variables: ***
  3849. ;;; mode: C++ ***
  3850. ;;; End: ***
  3851. */
  3852.